meshGFaceTransfinite.cpp

// Gmsh - Copyright (C) 1997-2013 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.

#include <map>
#include "meshGFace.h"
#include "GVertex.h"
#include "GEdge.h"
#include "GFace.h"
#include "MVertex.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "Context.h"
#include "GmshMessage.h"
#include "Numeric.h"

#define SQU(a)      ((a)*(a))

/*
   s4 +-----c3-----+ s3
      |            |
      |            |
     c4            c2
      |            |
      |            |
   s1 +-----c1-----+ s2
*/

// f(u,v) = (1-u) c4(v) + u c2(v) + (1-v) c1(u) + v c3(u)
//          - [ (1-u)(1-v) s1 + u(1-v) s2 + uv s3 + (1-u)v s4 ]
#define TRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
   (1.-u)*c4+u*c2+(1.-v)*c1+v*c3-((1.-u)*(1.-v)*s1+u*(1.-v)*s2+u*v*s3+(1.-u)*v*s4)

// s1=s4=c4
// f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3
#define TRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) u*c2+(1.-v)*c1+v*c3-(u*(1.-v)*s2+u*v*s3)

void findTransfiniteCorners(GFace *gf, std::vector<MVertex*> &corners)
{
  if(gf->meshAttributes.corners.size()){
    // corners have been specified explicitly
    for(unsigned int i = 0; i < gf->meshAttributes.corners.size(); i++)
      corners.push_back(gf->meshAttributes.corners[i]->mesh_vertices[0]);
  }
  else{
    // try to find the corners automatically
    std::list<GEdge*> fedges = gf->edges();
    GEdgeLoop el(fedges);
    for(GEdgeLoop::iter it = el.begin(); it != el.end(); it++)
      corners.push_back(it->getBeginVertex()->mesh_vertices[0]);

    // try reaaally hard for 3-sided faces
    if(corners.size() == 3){
      GEdge *first = 0, *last = 0;
      for(std::list<GEdge*>::iterator it = fedges.begin(); it != fedges.end(); it++){
        if(((*it)->getBeginVertex()->mesh_vertices[0] == corners[0] &&
            (*it)->getEndVertex()->mesh_vertices[0] == corners[1]) ||
           ((*it)->getBeginVertex()->mesh_vertices[0] == corners[1] &&
            (*it)->getEndVertex()->mesh_vertices[0] == corners[0])){
          first = *it;
        }
        if(((*it)->getBeginVertex()->mesh_vertices[0] == corners[2] &&
            (*it)->getEndVertex()->mesh_vertices[0] == corners[0]) ||
           ((*it)->getBeginVertex()->mesh_vertices[0] == corners[0] &&
            (*it)->getEndVertex()->mesh_vertices[0] == corners[2])){
          last = *it;
        }
      }
      if(first && last){