Skip to content
Snippets Groups Projects
Generator.cpp 37.2 KiB
Newer Older
Christophe Geuzaine's avatar
Christophe Geuzaine committed
// Gmsh - Copyright (C) 1997-2017 C. Geuzaine, J.-F. Remacle
Christophe Geuzaine's avatar
Christophe Geuzaine committed
//
Christophe Geuzaine's avatar
Christophe Geuzaine committed
// See the LICENSE.txt file for license information. Please report all
Christophe Geuzaine's avatar
Christophe Geuzaine committed
// bugs and problems to the public mailing list <gmsh@onelab.info>.
Claudine Bon's avatar
Claudine Bon committed
#include <stdlib.h>
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
#include "OS.h"
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
#include "GModel.h"
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
#include "MPoint.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
#include "meshGEdge.h"
#include "meshGFace.h"
#include "meshGFaceOptimize.h"
#include "meshGFaceBDS.h"
#include "meshGRegion.h"
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
#include "BackgroundMesh.h"
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
#include "HighOrder.h"
#include "Generator.h"
Christophe Geuzaine's avatar
Christophe Geuzaine committed
#include "meshGFaceLloyd.h"
#include "GFaceCompound.h"
#include "Options.h"
#include "simple3D.h"
#include "meshGRegionRelocateVertex.h"
#if defined(HAVE_OPTHOM)
#include "OptHomRun.h"
#include "OptHomElastic.h"
#endif

#if defined(HAVE_POST)
class TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES {
public:
  void operator () (GRegion *gr) {
    std::list<GEdge*> e = gr->embeddedEdges();
    std::list<GFace*> f = gr->embeddedFaces();
    if (e.empty() && f.empty())return;
    std::map<MEdge,GEdge*,Less_Edge> edges;
    std::map<MFace,GFace*,Less_Face> faces;
    std::list<GEdge*>::iterator it = e.begin();
    std::list<GFace*>::iterator itf = f.begin();
    for ( ; it != e.end() ; ++it){
      for (unsigned int i=0;i<(*it)->lines.size(); ++i){
	if (distance ((*it)->lines[i]->getVertex(0),(*it)->lines[i]->getVertex(1)) > 1.e-12)
	  edges.insert(std::make_pair(MEdge((*it)->lines[i]->getVertex(0),(*it)->lines[i]->getVertex(1)),*it));
      }
    }
    for ( ; itf != f.end() ; ++itf){
      for (unsigned int i=0;i<(*itf)->triangles.size(); ++i){
	faces.insert(std::make_pair(MFace((*itf)->triangles[i]->getVertex(0),(*itf)->triangles[i]->getVertex(1),(*itf)->triangles[i]->getVertex(2)),*itf));
      }
    }
    Msg::Info ("Searching for %d embedded mesh edges and %d embedded mesh faces in region %d", edges.size(),  faces.size(), gr->tag());
    for (unsigned int k=0;k<gr->getNumMeshElements();k++){
      for (int j=0;j<gr->getMeshElement(k)->getNumEdges();j++){
	edges.erase (gr->getMeshElement(k)->getEdge(j));
      }
      for (int j=0;j<gr->getMeshElement(k)->getNumFaces();j++){
	faces.erase (gr->getMeshElement(k)->getFace(j));
      }
    }
    if (edges.empty() && faces.empty()) {
      Msg::Info ("All embedded edges and faces are present in the final mesh");
    }
    if (edges.size()) {
      char name[256];
      sprintf(name,"missingEdgesOnRegion%d.pos",gr->tag());
      Msg::Error("Region %d : %d mesh edges that should be embedded are missing in the final mesh",gr->tag(), (int)edges.size());
      Msg::Error("Saving the missing edges in file %s",name);
      FILE *f = fopen(name,"w");
      fprintf(f,"View \" \" {\n");
      for (std::map<MEdge,GEdge*,Less_Edge>::iterator it =  edges.begin() ; it != edges.end(); ++it){
	MVertex *v1 = it->first.getVertex(0);
	MVertex *v2 = it->first.getVertex(1);
	fprintf(f,"SL(%g,%g,%g,%g,%g,%g){%d,%d};\n",v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(), it->second->tag(),it->second->tag());
      }
      fprintf(f,"};\n");
      fclose(f);
    }
    if (faces.size()) {
      char name[256];
      sprintf(name,"missingFacesOnRegion%d.pos",gr->tag());
      Msg::Error("Region %d : %d mesh faces that should be embedded are missing in the final mesh",gr->tag(), (int)faces.size());
      Msg::Error("Saving the missing faces in file %s",name);
      FILE *f = fopen(name,"w");
      fprintf(f,"View \" \" {\n");
      for (std::map<MFace,GFace*,Less_Face>::iterator it =  faces.begin() ; it != faces.end(); ++it){
	MVertex *v1 = it->first.getVertex(0);
	MVertex *v2 = it->first.getVertex(1);
	MVertex *v3 = it->first.getVertex(2);
	fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(),
		v3->x(),v3->y(),v3->z(),it->second->tag(),it->second->tag(),it->second->tag());
      }
      fprintf(f,"};\n");
      fclose(f);
    }
  }
};

static void GetQualityMeasure(std::vector<T*> &ele,
                              double &gamma, double &gammaMin, double &gammaMax,
                              double &minSICN, double &minSICNMin, double &minSICNMax,
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
                              double &rho, double &rhoMin, double &rhoMax,
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
{
  for(unsigned int i = 0; i < ele.size(); i++){
    double g = ele[i]->gammaShapeMeasure();
    gamma += g;
    gammaMin = std::min(gammaMin, g);
    gammaMax = std::max(gammaMax, g);
    double s = ele[i]->minSICNShapeMeasure();
    minSICN += s;
    minSICNMin = std::min(minSICNMin, s);
    minSICNMax = std::max(minSICNMax, s);
    double r = ele[i]->rhoShapeMeasure();
    rho += r;
    rhoMin = std::min(rhoMin, r);
    rhoMax = std::max(rhoMax, r);
    for(int j = 0; j < 100; j++){
      if(s > (2*j-100) / 100. && s <= (2*j-98) / 100.) quality[0][j]++;
      if(g > j / 100. && g <= (j + 1) / 100.) quality[1][j]++;
      if(r > j / 100. && r <= (j + 1) / 100.) quality[2][j]++;
    }
  }
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
}

void GetStatistics(double stat[50], double quality[3][100])
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
{
  for(int i = 0; i < 50; i++) stat[i] = 0.;

Christophe Geuzaine's avatar
Christophe Geuzaine committed
  GModel *m = GModel::current();

Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  if(!m) return;

  stat[0] = m->getNumVertices();
  stat[1] = m->getNumEdges();
  stat[2] = m->getNumFaces();
  stat[3] = m->getNumRegions();
  std::map<int, std::vector<GEntity*> > physicals[4];
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  m->getPhysicalGroups(physicals);
  stat[45] = physicals[0].size() + physicals[1].size() +
    physicals[2].size() + physicals[3].size();
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
    stat[4] += (*it)->mesh_vertices.size();
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
    // TODO: make this an option! if((*it)->getVisibility()){
    stat[5] += (*it)->mesh_vertices.size();
    stat[7] += (*it)->triangles.size();
    stat[8] += (*it)->quadrangles.size();
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
    stat[6] += (*it)->mesh_vertices.size();
    stat[9] += (*it)->tetrahedra.size();
    stat[10] += (*it)->hexahedra.size();
    stat[11] += (*it)->prisms.size();
    stat[12] += (*it)->pyramids.size();
    stat[13] += (*it)->trihedra.size();
  stat[14] = CTX::instance()->meshTimer[0];
  stat[15] = CTX::instance()->meshTimer[1];
  stat[16] = CTX::instance()->meshTimer[2];
  if(quality){
    for(int i = 0; i < 3; i++)
      for(int j = 0; j < 100; j++)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
        quality[i][j] = 0.;
    double minSICN = 0., minSICNMin = 1., minSICNMax = -1.;
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    double gamma = 0., gammaMin = 1., gammaMax = 0.;
    double rho = 0., rhoMin = 1., rhoMax = 0.;

    double N = stat[9] + stat[10] + stat[11] + stat[12] + stat[13];
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    if(N){ // if we have 3D elements
      for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
        GetQualityMeasure((*it)->tetrahedra, gamma, gammaMin, gammaMax,
                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
        GetQualityMeasure((*it)->hexahedra, gamma, gammaMin, gammaMax,
                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
        GetQualityMeasure((*it)->prisms, gamma, gammaMin, gammaMax,
                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
        GetQualityMeasure((*it)->pyramids, gamma, gammaMin, gammaMax,
                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
    }
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    else{ // 2D elements
      N = stat[7] + stat[8];
      for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
        GetQualityMeasure((*it)->quadrangles, gamma, gammaMin, gammaMax,
                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
        GetQualityMeasure((*it)->triangles, gamma, gammaMin, gammaMax,
                          minSICN, minSICNMin, minSICNMax, rho, rhoMin, rhoMax, quality);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
      }
    }
    if(N){
      stat[18] = minSICN / N; stat[19] = minSICNMin; stat[20] = minSICNMax;
      stat[21] = gamma / N;   stat[22] = gammaMin;   stat[23] = gammaMax;
      stat[25] = rho / N;     stat[25] = rhoMin;     stat[26] = rhoMax;
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    }
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

#if defined(HAVE_POST)
  stat[27] = PView::list.size();
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  for(unsigned int i = 0; i < PView::list.size(); i++) {
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    PViewData *data = PView::list[i]->getData(true);
    stat[28] += data->getNumPoints();
    stat[29] += data->getNumLines();
    stat[30] += data->getNumTriangles();
    stat[31] += data->getNumQuadrangles();
    stat[32] += data->getNumTetrahedra();
    stat[33] += data->getNumHexahedra();
    stat[34] += data->getNumPrisms();
    stat[35] += data->getNumPyramids();
    stat[36] += data->getNumStrings2D() + data->getNumStrings3D();
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  }
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
static bool TooManyElements(GModel *m, int dim)
  if(CTX::instance()->expertMode || !m->getNumVertices()) return false;

  // try to detect obvious mistakes in characteristic lenghts (one of
  // the most common cause for erroneous bug reports on the mailing
  // list)
  double sumAllLc = 0.;
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it)
    sumAllLc += (*it)->prescribedMeshSizeAtVertex() * CTX::instance()->mesh.lcFactor;
  sumAllLc /= (double)m->getNumVertices();
  if(!sumAllLc || pow(CTX::instance()->lc / sumAllLc, dim) > 1.e10)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
    return !Msg::GetAnswer
      ("Your choice of mesh element sizes will likely produce a very\n"
       "large mesh. Do you really want to continue?\n\n"
       "(To disable this warning in the future, select `Enable expert mode'\n"
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
       "in the option dialog.)", 1, "Cancel", "Continue");
  return false;
}

static bool CancelDelaunayHybrid(GModel *m)
{
  if(CTX::instance()->expertMode) return false;
  int n = 0;
  for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
    n += (*it)->getNumMeshElements();
  if(n)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
    return !Msg::GetAnswer
      ("You are trying to generate a mixed structured/unstructured grid using\n"
       "the 3D Delaunay algorithm. This algorithm cannot garantee that the\n"
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
       "final mesh will be conforming. (You should probably use the 3D Frontal\n"
       "algorithm instead.) Do you really want to continue with the Delaunay?\n\n"
       "(To disable this warning in the future, select `Enable expert mode'\n"
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
       "in the option dialog.)", 1, "Cancel", "Continue");
  return false;
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
static void Mesh0D(GModel *m)
{
  for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it){
    GVertex *gv = *it;
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
    if(gv->mesh_vertices.empty())
      gv->mesh_vertices.push_back(new MVertex(gv->x(), gv->y(), gv->z(), gv));
    if(gv->points.empty())
      gv->points.push_back(new MPoint(gv->mesh_vertices.back()));
  }
  for(GModel::viter it = m->firstVertex(); it != m->lastVertex(); ++it){
    GVertex *gv = *it;
      if (gv->correspondingVertices.empty()){
        GVertex *master = dynamic_cast<GVertex*>(gv->meshMaster());
        if(master)
          gv->correspondingVertices[gv->mesh_vertices[0]] = master->mesh_vertices[0];
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
}

Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
static void Mesh1D(GModel *m)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
{
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  if(TooManyElements(m, 1)) return;
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

Christophe Geuzaine's avatar
Christophe Geuzaine committed
  for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
    (*it)->meshStatistics.status = GEdge::PENDING;

  Msg::ResetProgressMeter();

  int nIter = 0, nTot = m->getNumEdges();
    int nPending = 0;
    for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it){
Christophe Geuzaine's avatar
Christophe Geuzaine committed
      if ((*it)->meshStatistics.status == GEdge::PENDING){
      if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 1D...");
    if(!nPending) break;
    if(nIter++ > 10) break;
  }

  CTX::instance()->meshTimer[0] = t2 - t1;
  Msg::StatusBar(true, "Done meshing 1D (%g s)", CTX::instance()->meshTimer[0]);
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
static void PrintMesh2dStatistics(GModel *m)
  FILE *statreport = 0;
  if(CTX::instance()->createAppendMeshStatReport == 1)
    statreport = Fopen(CTX::instance()->meshStatReportFileName.c_str(), "w");
  else if(CTX::instance()->createAppendMeshStatReport == 2)
    statreport = Fopen(CTX::instance()->meshStatReportFileName.c_str(), "a");
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  else
    return;

Christophe Geuzaine's avatar
Christophe Geuzaine committed
  if(!statreport){
    Msg::Error("Could not open file '%s'",
               CTX::instance()->meshStatReportFileName.c_str());
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    return;
  }
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  double worst = 1, best = 0, avg = 0;
  double e_long = 0, e_short = 1.e22, e_avg = 0;
  int nTotT = 0, nTotE = 0, nTotGoodLength = 0, nTotGoodQuality = 0;
  int nUnmeshed = 0, numFaces = 0;
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  if(CTX::instance()->createAppendMeshStatReport == 1){
    fprintf(statreport, "2D stats\tname\t\t#faces\t\t#fail\t\t"
            "#t\t\tQavg\t\tQbest\t\tQworst\t\t#Q>90\t\t#Q>90/#t\t"
            "#e\t\ttau\t\t#Egood\t\t#Egood/#e\tCPU\n");
    if(m->empty()){
      fclose(statreport);
      return;
    }
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  }

Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  for(GModel::fiter it = m->firstFace() ; it != m->lastFace(); ++it){
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
    if((*it)->geomType() != GEntity::DiscreteSurface){
      worst = std::min((*it)->meshStatistics.worst_element_shape, worst);
      best = std::max((*it)->meshStatistics.best_element_shape, best);
      avg += (*it)->meshStatistics.average_element_shape *
        (*it)->meshStatistics.nbTriangle;
      e_avg += (*it)->meshStatistics.efficiency_index;
      e_long = std::max((*it)->meshStatistics.longest_edge_length, e_long);
      e_short = std::min((*it)->meshStatistics.smallest_edge_length, e_short);
      if ((*it)->meshStatistics.status == GFace::FAILED ||
          (*it)->meshStatistics.status == GFace::PENDING) nUnmeshed++;
      nTotT += (*it)->meshStatistics.nbTriangle;
      nTotE += (*it)->meshStatistics.nbEdge;
      nTotGoodLength += (*it)->meshStatistics.nbGoodLength;
      nTotGoodQuality += (*it)->meshStatistics.nbGoodQuality;
      numFaces++;
    }
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  }
  Msg::Info("*** Efficiency index for surface mesh tau=%g ",
            100*exp(e_avg/(double)nTotE));
  fprintf(statreport,"\t%16s\t%d\t\t%d\t\t", m->getName().c_str(), numFaces,
          nUnmeshed);
  fprintf(statreport,"%d\t\t%8.7f\t%8.7f\t%8.7f\t%d\t\t%8.7f\t",
          nTotT, avg / (double)nTotT, best, worst, nTotGoodQuality,
          (double)nTotGoodQuality / nTotT);
  fprintf(statreport,"%d\t\t%8.7f\t%d\t\t%8.7f\t%8.1f\n",
          nTotE, exp(e_avg / (double)nTotE), nTotGoodLength,
          (double)nTotGoodLength / nTotE, CTX::instance()->meshTimer[1]);
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  fclose(statreport);
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
static void Mesh2D(GModel *m)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
{
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  if(TooManyElements(m, 2)) return;
  double t1 = GetTimeInSeconds();
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

Christophe Geuzaine's avatar
Christophe Geuzaine committed
  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
    (*it)->meshStatistics.status = GFace::PENDING;

  // boundary layers are special: their generation (including vertices and curve
  // meshes) is global as it depends on a smooth normal field generated from the
  // surface mesh of the source surfaces
  if(!Mesh2DWithBoundaryLayers(m)){
    for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
      if ((*it)->geomType() == GEntity::CompoundSurface)
        cf.insert(*it);
      else
        f.insert(*it);
    Msg::ResetProgressMeter();

    int nIter = 0, nTot = m->getNumFaces();
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
    while(1){
      int nPending = 0;
      std::vector<GFace*> temp;
      temp.insert(temp.begin(), f.begin(), f.end());
Gaetan Bricteux's avatar
Gaetan Bricteux committed
#if defined(_OPENMP)
#pragma omp parallel for schedule (dynamic)
Gaetan Bricteux's avatar
Gaetan Bricteux committed
#endif
      for(size_t K = 0 ; K < temp.size() ; K++){
        if (temp[K]->meshStatistics.status == GFace::PENDING){
Christophe Geuzaine's avatar
Christophe Geuzaine committed
          backgroundMesh::current()->unset();
Christophe Geuzaine's avatar
Christophe Geuzaine committed
          if(CTX::instance()->mesh.optimizeLloyd){
            if (temp[K]->geomType()==GEntity::CompoundSurface ||
                temp[K]->geomType()==GEntity::Plane ||
                temp[K]->geomType()==GEntity::RuledSurface) {
Christophe Geuzaine's avatar
Christophe Geuzaine committed
              if (temp[K]->meshAttributes.method != MESH_TRANSFINITE &&
                  !temp[K]->meshAttributes.extrude) {
                smoothing smm(CTX::instance()->mesh.optimizeLloyd, 6);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
                //m->writeMSH("beforeLLoyd.msh");
                smm.optimize_face(temp[K]);
                int rec = ((CTX::instance()->mesh.recombineAll ||
                            temp[K]->meshAttributes.recombine) &&
                           !CTX::instance()->mesh.recombine3DAll);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
                //m->writeMSH("afterLLoyd.msh");
                if (rec) recombineIntoQuads(temp[K]);
                //m->writeMSH("afterRecombine.msh");
              }
            }
          }
Gaetan Bricteux's avatar
Gaetan Bricteux committed
#if defined(_OPENMP)
#pragma omp critical
Gaetan Bricteux's avatar
Gaetan Bricteux committed
#endif
        if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D...");
Gaetan Bricteux's avatar
Gaetan Bricteux committed
#if defined(_OPENMP)
#pragma omp master
Gaetan Bricteux's avatar
Gaetan Bricteux committed
#endif
      for(std::set<GFace*, GEntityLessThan>::iterator it = cf.begin();
          it != cf.end(); ++it){
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
        if ((*it)->meshStatistics.status == GFace::PENDING){
Christophe Geuzaine's avatar
Christophe Geuzaine committed
          backgroundMesh::current()->unset();
	  //          meshGFace mesher(true);
          (*it)->mesh(true);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
          if(CTX::instance()->mesh.optimizeLloyd){
            if ((*it)->geomType()==GEntity::CompoundSurface ||
                (*it)->geomType()==GEntity::Plane ||
                (*it)->geomType()==GEntity::RuledSurface) {
Christophe Geuzaine's avatar
Christophe Geuzaine committed
              if ((*it)->meshAttributes.method != MESH_TRANSFINITE &&
                  !(*it)->meshAttributes.extrude) {
                smoothing smm(CTX::instance()->mesh.optimizeLloyd, 6);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
                //m->writeMSH("beforeLLoyd.msh");
                smm.optimize_face(*it);
                int rec = ((CTX::instance()->mesh.recombineAll ||
                            (*it)->meshAttributes.recombine) &&
                           !CTX::instance()->mesh.recombine3DAll);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
                //m->writeMSH("afterLLoyd.msh");
                if (rec) recombineIntoQuads(*it);
                //m->writeMSH("afterRecombine.msh");
              }
            }
          }
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
        }
        if(!nIter) Msg::ProgressMeter(nPending, nTot, false, "Meshing 2D...");
      if(!nPending) break;
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
      if(nIter++ > 10) break;
  // collapseSmallEdges(*m);
  double t2 = GetTimeInSeconds();
  CTX::instance()->meshTimer[1] = t2 - t1;
  Msg::StatusBar(true, "Done meshing 2D (%g s)", CTX::instance()->meshTimer[1]);

  PrintMesh2dStatistics(m);
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
static void FindConnectedRegions(const std::vector<GRegion*> &del,
                                 std::vector<std::vector<GRegion*> > &connected)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
{
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
  std::vector<GRegion*> delaunay = del;
  // test: connected.resize(1); connected[0] = delaunay; return;

  const unsigned int nbVolumes = delaunay.size();
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
  if (!nbVolumes) return;
  while (delaunay.size()){
    std::set<GRegion*> oneDomain;
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
    std::stack<GRegion*> _stack;
    GRegion *r = delaunay[0];
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
    _stack.push(r);
    while(!_stack.empty()){
      r = _stack.top();
      _stack.pop();
      oneDomain.insert(r);
      std::list<GFace*> faces = r->faces();
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
      for (std::list<GFace*>::iterator it = faces.begin(); it != faces.end() ; ++it){
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
        GRegion *other = (gf->getRegion(0) == r) ? gf->getRegion(1) : gf->getRegion(0);
        if(other != 0 && oneDomain.find(other) == oneDomain.end())
          _stack.push(other);
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
    std::vector<GRegion*> temp1, temp2;
    for (unsigned int i = 0; i < delaunay.size(); i++){
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
      if (oneDomain.find(r) == oneDomain.end()) temp1.push_back(r);
      else temp2.push_back(r);
    }
    connected.push_back(temp2);
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
    delaunay = temp1;
  }
  Msg::Info("Delaunay Meshing %d volumes with %d connected components",
Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
            nbVolumes, connected.size());
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
}

template <class ITERATOR>
void fillv_(std::multimap<MVertex*, MElement*> &vertexToElement,
	    ITERATOR it_beg, ITERATOR it_end)
{
  for (ITERATOR IT = it_beg; IT != it_end ; ++IT){
    MElement *el = *IT;
    for(int j = 0; j < el->getNumVertices(); j++) {
      MVertex* e = el->getVertex(j);
      vertexToElement.insert(std::make_pair(e, el));
    }
  }
}

int LaplaceSmoothing(GRegion *gr)
  std::multimap<MVertex*, MElement*> vertexToElement;
  fillv_(vertexToElement, (gr)->tetrahedra.begin(), (gr)->tetrahedra.end());
  fillv_(vertexToElement, (gr)->hexahedra.begin(),  (gr)->hexahedra.end());
  fillv_(vertexToElement, (gr)->prisms.begin(),     (gr)->prisms.end());
  fillv_(vertexToElement, (gr)->pyramids.begin(),   (gr)->pyramids.end());
  int N=0;
  for (unsigned int i=0; i<gr->mesh_vertices.size();i++){
    MVertex *v = gr->mesh_vertices[i];
    std::multimap<MVertex*, MElement*>::iterator it = vertexToElement.lower_bound(v);
    std::multimap<MVertex*, MElement*>::iterator it_low = it;
    std::multimap<MVertex*, MElement*>::iterator it_up  = vertexToElement.upper_bound(v);
    double minQual = 1.e22;
    double volTot = 0.0;
    double xold=v->x(), yold=v->y(), zold=v->z();
    SPoint3 pNew (0,0,0);
    for (; it != it_up; ++it) {
      minQual= std::min(minQual,it->second->minSICNShapeMeasure());
      double vol = fabs(it->second->getVolume());
      SPoint3 cog = it->second->barycenter();
      pNew += cog * vol;
      volTot += vol;
    }
    pNew *= (1./volTot);
    v->setXYZ (pNew.x(),pNew.y(),pNew.z());
    double minQual2 = 1.e22;
    for (it = it_low; it != it_up; ++it) {
      minQual2 = std::min(minQual2,it->second->minSICNShapeMeasure());
	v->setXYZ (xold,yold,zold);
	break;
      }
    }
    if (minQual < minQual2) N++;
  return N;
}

// JFR : use hex-splitting to resolve non conformity
//     : if howto == 1 ---> split hexes
//     : if howto == 2 ---> create transition elements

/*
 v3        v2
  x--------x
  |\       |
  |  \     |
  |    \   |
  |      \ |
  x--------x
  v0       v1
 */

void buildUniqueFaces (GRegion *gr, std::set<MFace,Less_Face> &bnd)
{
  for (unsigned int i=0;i<gr->getNumMeshElements();i++){
    MElement *e = gr->getMeshElement(i);
    for(int j=0;j<e->getNumFaces();j++){
      MFace f = e->getFace(j);
      std::set<MFace,Less_Face>::iterator it = bnd.find(f);
      if (it == bnd.end())bnd.insert(f);
      else bnd.erase(it);
    }
  }
}

bool MakeMeshConformal(GModel *gm, int howto)
{
  fs_cont search;
  buildFaceSearchStructure(gm, search);
  std::set<MFace,Less_Face> bnd;
  for (GModel::riter rit = gm->firstRegion(); rit != gm->lastRegion(); ++rit){
    GRegion *gr = *rit;
    buildUniqueFaces (gr,bnd);
  }
  // bnd2 contains non conforming faces
   std::set<MFace,Less_Face> bnd2;
  for (std::set<MFace,Less_Face>::iterator itf = bnd.begin(); itf != bnd.end(); ++itf){
    GFace *gfound = findInFaceSearchStructure (*itf,search);
    if (!gfound){
      bnd2.insert(*itf);
    }
  }
  bnd.clear();

  Msg::Info("%d hanging faces",bnd2.size());

  std::set<MFace,Less_Face> ncf;
  for (std::set<MFace,Less_Face>::iterator itf = bnd2.begin(); itf != bnd2.end(); ++itf){
    const MFace &f = *itf;
    if (f.getNumVertices() == 4){ // quad face
      std::set<MFace,Less_Face>::iterator it1 =
        bnd2.find (MFace(f.getVertex(0),f.getVertex(1),f.getVertex(2)));
      std::set<MFace,Less_Face>::iterator it2 =
        bnd2.find (MFace(f.getVertex(2),f.getVertex(3),f.getVertex(0)));
      if (it1 != bnd2.end() && it2 != bnd2.end()){
	ncf.insert(MFace(f.getVertex(1), f.getVertex(2),
                         f.getVertex(3),f.getVertex(0)));
      }
      else {
	it1 = bnd2.find (MFace(f.getVertex(0),f.getVertex(1),f.getVertex(3)));
	it2 = bnd2.find (MFace(f.getVertex(3),f.getVertex(1),f.getVertex(2)));
	if (it1 != bnd2.end() && it2 != bnd2.end()){
	  ncf.insert(MFace(f.getVertex(0), f.getVertex(1),
                           f.getVertex(2), f.getVertex(3)));
	  Msg::Error("MakeMeshConformal: wrong mesh topology");
	  return false;
	}
      }
    }
  }
  bnd2.clear();

  for (GModel::riter rit = gm->firstRegion(); rit != gm->lastRegion(); ++rit){
    GRegion *gr = *rit;
    std::vector<MHexahedron*> remainingHexes;
    for (unsigned int i=0;i<gr->hexahedra.size();i++){
      MHexahedron *e = gr->hexahedra[i];
      std::vector<MFace> faces;
      for(int j=0;j<e->getNumFaces();j++){
	MFace f = e->getFace(j);
	std::set<MFace,Less_Face>::iterator it = ncf.find(f);
	if (it == ncf.end()){
	  faces.push_back(f);
	}
	else {
	  faces.push_back(MFace(it->getVertex(0),it->getVertex(1),it->getVertex(3)));
	  faces.push_back(MFace(it->getVertex(1),it->getVertex(2),it->getVertex(3)));
	}
      }
      // HEX IS ONLY SURROUNED BY COMPATIBLE ELEMENTS
      if ((int)faces.size() == e->getNumFaces()){
	remainingHexes.push_back(e);
      }
      else {
	SPoint3 pp = e->barycenter();
	MVertex *newv = new MVertex (pp.x(),pp.y(),pp.z(),gr);
	gr->mesh_vertices.push_back(newv);
	for (unsigned int j=0;j<faces.size();j++){
	  MFace &f = faces[j];
	  if (f.getNumVertices() == 4){
	    gr->pyramids.push_back(new MPyramid(f.getVertex(0), f.getVertex(1),
                                                f.getVertex(2), f.getVertex(3), newv));
	    gr->tetrahedra.push_back(new MTetrahedron(f.getVertex(0), f.getVertex(1),
                                                      f.getVertex(2), newv));
    }
    gr->hexahedra = remainingHexes;
    remainingHexes.clear();
    std::vector<MPrism*> remainingPrisms;
    for (unsigned int i=0;i<gr->prisms.size();i++){
      MPrism *e = gr->prisms[i];
      std::vector<MFace> faces;
      for(int j=0;j<e->getNumFaces();j++){
	MFace f = e->getFace(j);
	std::set<MFace,Less_Face>::iterator it = ncf.find(f);
	if (it == ncf.end()){
	  faces.push_back(f);
	}
	else {
	  faces.push_back(MFace(it->getVertex(0),it->getVertex(1),it->getVertex(3)));
	  faces.push_back(MFace(it->getVertex(1),it->getVertex(2),it->getVertex(3)));
	}
      }
      // HEX IS ONLY SURROUNED BY COMPATIBLE ELEMENTS
      if ((int)faces.size() == e->getNumFaces()){
	remainingPrisms.push_back(e);
      }
      else {
	SPoint3 pp = e->barycenter();
	MVertex *newv = new MVertex (pp.x(),pp.y(),pp.z(),gr);
	gr->mesh_vertices.push_back(newv);
	for (unsigned int j=0;j<faces.size();j++){
	  MFace &f = faces[j];
	  if (f.getNumVertices() == 4){
	    gr->pyramids.push_back(new MPyramid(f.getVertex(0), f.getVertex(1),
                                                f.getVertex(2), f.getVertex(3), newv));
	    gr->tetrahedra.push_back(new MTetrahedron(f.getVertex(0), f.getVertex(1),
                                                      f.getVertex(2), newv));
    }
    gr->prisms = remainingPrisms;
  }

  return true;
}

  fs_cont search;
  buildFaceSearchStructure(gm, search);
  int count = 0;
  for (GModel::riter rit = gm->firstRegion(); rit != gm->lastRegion(); ++rit){
    GRegion *gr = *rit;
    std::set<MFace,Less_Face> bnd;
    double vol = 0.0;
    for (unsigned int i=0;i<gr->getNumMeshElements();i++){
      MElement *e = gr->getMeshElement(i);
      vol += fabs(e->getVolume());
      for(int j=0;j<e->getNumFaces();j++){
	MFace f = e->getFace(j);
	std::set<MFace,Less_Face>::iterator it = bnd.find(f);
	if (it == bnd.end())bnd.insert(f);
	else bnd.erase(it);
      }
    }
    printf("vol(%d) = %12.5E\n",gr->tag(),vol);
    for (std::set<MFace,Less_Face>::iterator itf = bnd.begin(); itf != bnd.end(); ++itf){
      GFace *gfound = findInFaceSearchStructure (*itf,search);
      if (!gfound){
	count ++;
      }
    }
  }
  if (!count)Msg::Info("Mesh Conformity: OK");
  else Msg::Error ("Mesh is not conforming (%d hanging faces)!",count);
}

Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
static void Mesh3D(GModel *m)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
{
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  if(TooManyElements(m, 3)) return;
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

  // mesh the extruded volumes first
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  std::for_each(m->firstRegion(), m->lastRegion(), meshGRegionExtruded());
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  // then subdivide if necessary (unfortunately the subdivision is a
  // global operation, which can require changing the surface mesh!)
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  SubdivideExtrudedMesh(m);
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

  // then mesh all the non-delaunay regions (front3D with netgen)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  std::vector<GRegion*> delaunay;
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  std::for_each(m->firstRegion(), m->lastRegion(), meshGRegion(delaunay));
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

  // warn if attempting to use Delaunay for mixed meshes
  if(delaunay.size() && CancelDelaunayHybrid(m)) return;

Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
  // and finally mesh the delaunay regions (again, this is global; but
  // we mesh each connected part separately for performance and mesh
  // quality reasons)
  std::vector<std::vector<GRegion*> > connected;
  FindConnectedRegions(delaunay, connected);
  // remove quads elements for volumes that are recombined
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  // pragma OMP here ?
  for(unsigned int i = 0; i < connected.size(); i++){
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    for(unsigned j = 0; j < connected[i].size(); j++){
      GRegion *gr = connected[i][j];
      if(CTX::instance()->mesh.recombine3DAll || gr->meshAttributes.recombine3D){
Christophe Geuzaine's avatar
Christophe Geuzaine committed
        for(std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it)
          quadsToTriangles(*it, 1000000);
Christophe Geuzaine's avatar
Christophe Geuzaine committed

Christophe Geuzaine's avatar
Christophe Geuzaine committed
  double time_recombination = 0., vol_element_recombination = 0.;
  double vol_hexa_recombination = 0.;
  int nb_elements_recombination = 0, nb_hexa_recombination = 0;

  // pragma OMP here ?
Christophe Geuzaine's avatar
Christophe Geuzaine committed
  for(unsigned int i = 0; i < connected.size(); i++){
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
    MeshDelaunayVolume(connected[i]);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    // additional code for experimental hex mesh
    for(unsigned j = 0; j < connected[i].size(); j++){
Emilie Marchandise's avatar
Emilie Marchandise committed
      GRegion *gr = connected[i][j];
Christophe Geuzaine's avatar
Christophe Geuzaine committed
      bool treat_region_ok = false;
Emilie Marchandise's avatar
Emilie Marchandise committed
      if(CTX::instance()->mesh.algo3d == ALGO_3D_RTREE){
Christophe Geuzaine's avatar
Christophe Geuzaine committed
          treat_region_ok = true;
Emilie Marchandise's avatar
Emilie Marchandise committed
      }
Christophe Geuzaine's avatar
Christophe Geuzaine committed
      if(treat_region_ok && (CTX::instance()->mesh.recombine3DAll ||
                             gr->meshAttributes.recombine3D)){
        if (CTX::instance()->mesh.optimize){
          optimizeMeshGRegionGmsh opt;
          opt(gr);
        }
Christophe Geuzaine's avatar
Christophe Geuzaine committed
        double a = Cpu();
	//	CTX::instance()->mesh.recombine3DLevel = 2;
        if (CTX::instance()->mesh.recombine3DLevel >= 0){
          Recombinator rec;
          rec.execute(gr);
        }
        if (CTX::instance()->mesh.recombine3DLevel >= 1){
          Supplementary sup;
          sup.execute(gr);
        }
	printf("-----------> %d %d\n",CTX::instance()->mesh.recombine3DLevel,CTX::instance()->mesh.recombine3DConformity);
	post.execute(gr,CTX::instance()->mesh.recombine3DLevel,CTX::instance()->mesh.recombine3DConformity);
	//			     CTX::instance()->mesh.recombine3DConformity);
        // 0: no pyramid, 1: single-step, 2: two-steps (conforming),
        // true: fill non-conformities with trihedra
	RelocateVertices(gr, CTX::instance()->mesh.nbSmoothing);
        // while(LaplaceSmoothing (gr)){
        // }
	nb_elements_recombination += post.get_nb_elements();
	nb_hexa_recombination += post.get_nb_hexahedra();
	vol_element_recombination += post.get_vol_elements();
	vol_hexa_recombination += post.get_vol_hexahedra();
Christophe Geuzaine's avatar
Christophe Geuzaine committed
        // partial export
        //        stringstream ss;
        //        ss << "yamakawa_part_";
        //        ss << gr->tag();
        //        ss << ".msh";
        //        export_gregion_mesh(gr, ss.str().c_str());
Christophe Geuzaine's avatar
Christophe Geuzaine committed
        time_recombination += (Cpu() - a);
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    }
Christophe Geuzaine's avatar
Christophe Geuzaine committed
    Msg::Info("RECOMBINATION timing:");
    Msg::Info(" --- CUMULATIVE TIME RECOMBINATION : %g s.", time_recombination);
    Msg::Info("RECOMBINATION CUMULATIVE STATISTICS:");
    Msg::Info(".... Percentage of hexahedra   (#) : %g",
              nb_hexa_recombination*100./nb_elements_recombination);
    Msg::Info(".... Percentage of hexahedra (Vol) : %g",
              vol_hexa_recombination*100./vol_element_recombination);
    //    MakeMeshConformal (m, 1);
    TestConformity(m);
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

  // ensure that all volume Jacobians are positive
  //  std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionNetgen());
    std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
  CTX::instance()->mesh.changed = ENT_ALL;
  CTX::instance()->meshTimer[2] = t2 - t1;
  Msg::StatusBar(true, "Done meshing 3D (%g s)", CTX::instance()->meshTimer[2]);
void OptimizeMeshNetgen(GModel *m)
{
  Msg::StatusBar(true, "Optimizing 3D mesh with Netgen...");
  double t1 = Cpu();

  std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionNetgen());
  // Ensure that all volume Jacobians are positive
  m->setAllVolumesPositive();

    std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
  double t2 = Cpu();
  Msg::StatusBar(true, "Done optimizing 3D mesh with Netgen (%g s)", t2 - t1);
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
void OptimizeMesh(GModel *m)
  Msg::StatusBar(true, "Optimizing 3D mesh...");
  double t1 = Cpu();
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed

  std::for_each(m->firstRegion(), m->lastRegion(), optimizeMeshGRegionGmsh());
  // Ensure that all volume Jacobians are positive
  m->setAllVolumesPositive();

    std::for_each(m->firstRegion(), m->lastRegion(), TEST_IF_MESH_IS_COMPATIBLE_WITH_EMBEDDED_ENTITIES ());
  CTX::instance()->mesh.changed = ENT_ALL;
  double t2 = Cpu();
  Msg::StatusBar(true, "Done optimizing 3D mesh (%g s)", t2 - t1);
void SmoothMesh(GModel *m)
{
  Msg::StatusBar(true, "Smoothing 2D mesh...");
  double t1 = Cpu();

  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
    GFace *gf = *it;
    laplaceSmoothing(gf);
  }

  CTX::instance()->mesh.changed = ENT_ALL;
  double t2 = Cpu();
  Msg::StatusBar(true, "Done smoothing 2D mesh (%g s)", t2 - t1);
}

Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
void AdaptMesh(GModel *m)
  Msg::StatusBar(true, "Adapting 3D mesh...");
  double t1 = Cpu();

Christophe Geuzaine's avatar
pp  
Christophe Geuzaine committed
  for(int i = 0; i < 10; i++)
    std::for_each(m->firstRegion(), m->lastRegion(), adaptMeshGRegion());
  CTX::instance()->mesh.changed = ENT_ALL;
  double t2 = Cpu();
  Msg::StatusBar(true, "Done adaptating 3D mesh (%g s)", t2 - t1);
void RecombineMesh(GModel *m)
{
  Msg::StatusBar(true, "Recombining 2D mesh...");
  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
    GFace *gf = *it;
    recombineIntoQuads(gf);
  }

  CTX::instance()->mesh.changed = ENT_ALL;
  double t2 = Cpu();
  Msg::StatusBar(true, "Done recombining 2D mesh (%g s)", t2 - t1);
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
void GenerateMesh(GModel *m, int ask)
Christophe Geuzaine's avatar
 
Christophe Geuzaine committed
{
  // ProfilerStart("gmsh.prof");
  if(CTX::instance()->lock) {
    Msg::Info("I'm busy! Ask me that later...");