Forked from
gmsh / gmsh
8541 commits behind the upstream repository.
-
Christophe Geuzaine authoredChristophe Geuzaine authored
meshGFaceElliptic.cpp 27.55 KiB
// Gmsh - Copyright (C) 1997-2016 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@onelab.info>.
#include <stack>
#include "GmshConfig.h"
#include "meshGFaceElliptic.h"
#include "qualityMeasures.h"
#include "GFace.h"
#include "GEdge.h"
#include "GEdgeCompound.h"
#include "GVertex.h"
#include "GModel.h"
#include "MVertex.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MLine.h"
#include "BackgroundMeshTools.h"
#include "Numeric.h"
#include "GmshMessage.h"
#include "Generator.h"
#include "Context.h"
#include "OS.h"
#include "SVector3.h"
#include "SPoint3.h"
#include "fullMatrix.h"
#include "CenterlineField.h"
#if defined(HAVE_ANN)
#include "ANN/ANN.h"
#endif
#define TRAN_QUAD(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)
static void printQuads(GFace *gf, fullMatrix<SPoint2> uv,
std::vector<MQuadrangle*> quads, int iter)
{
if(!CTX::instance()->mesh.saveAll) return;
char name[234];
sprintf(name, "quadUV_%d_%d.pos", gf->tag(), iter);
FILE *f = Fopen(name, "w");
if(!f){
Msg::Error("Could not open file '%s'", name);
return;
}
fprintf(f,"View \"%s\" {\n",name);
for (int i = 1; i < uv.size1()-1; i++)
for (int j = 0; j < uv.size2(); j++)
fprintf(f,"SP(%g,%g,%g) {%d};\n", uv(i,j).x(), uv(i,j).y(), 0.0, i);
fprintf(f,"};\n");
fclose(f);
char name3[234];
sprintf(name3, "quadXYZ_%d_%d.pos", gf->tag(), iter);
FILE *f3 = Fopen(name3,"w");
if(!f3){
Msg::Error("Could not open file '%s'", name3);
return;
}
fprintf(f3,"View \"%s\" {\n",name3);
for (unsigned int i = 0; i < quads.size(); i++){
quads[i]->writePOS(f3,true,false,false,false,false,false);
}
fprintf(f3,"};\n"); fclose(f3);
}
static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<MVertex*> vert2,
std::vector<MVertex*> e01, std::vector<MVertex*> e10,
std::vector<MVertex*> e23, std::vector<MVertex*> e32,
std::vector<MVertex*> e02, std::vector<MVertex*> e13,
std::vector<MQuadrangle*> quads)
{
if(!CTX::instance()->mesh.saveAll) return;
std::vector<SPoint2> p1,p2;
for (unsigned int i = 0; i< vert1.size(); i++){
SPoint2 pi;
reparamMeshVertexOnFace(vert1[i], gf, pi);
p1.push_back(pi);
}
for (unsigned int j = 0; j< vert2.size(); j++){
SPoint2 pj;
reparamMeshVertexOnFace(vert2[j], gf, pj);
p2.push_back(pj);
}
char name[234];
sprintf(name,"paramGrid_%d.pos", gf->tag());
FILE *f = Fopen(name,"w");
if(!f){
Msg::Error("Could not open file '%s'", name);
return;
}
fprintf(f,"View \"%s\" {\n",name);
// for (unsigned int i = 0; i < p1.size(); i++)
// fprintf(f,"SP(%g,%g,%g) {%d};\n", p1[i].x(), p1[i].y(), 0.0, i);
// for (unsigned int j = 0; j < p2.size(); j++)
// fprintf(f,"SP(%g,%g,%g) {%d};\n", p2[j].x(), p2[j].y(), 0.0, 100+j);
for (unsigned int i = 0; i < p1.size()-1; i++)
fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p1[i].x(), p1[i].y(), 0.0, p1[i+1].x(), p1[i+1].y(), 0.0, 1, 1);
for (unsigned int i = 0; i < p2.size()-1; i++)
fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p2[i].x(), p2[i].y(), 0.0, p2[i+1].x(), p2[i+1].y(), 0.0, 1, 1);
fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p1[p1.size()-1].x(), p1[ p1.size()-1].y(), 0.0, p1[0].x(), p1[0].y(), 0.0, 1, 1);
fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p2[p2.size()-1].x(), p2[ p2.size()-1].y(), 0.0, p2[0].x(), p2[0].y(), 0.0, 1, 1);
fprintf(f,"};\n");
fclose(f);
char name2[234];
sprintf(name2,"paramEdges_%d.pos", gf->tag());
FILE *f2 = Fopen(name2,"w");
if(!f2){
Msg::Error("Could not open file '%s'", name2);
return;
}
fprintf(f2,"View \"%s\" {\n",name2);
for (unsigned int i = 0; i < e01.size(); i++){
SPoint2 pi; reparamMeshVertexOnFace(e01[i], gf, pi);
fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 1);
}
for (unsigned int i = 0; i < e10.size(); i++){
SPoint2 pi; reparamMeshVertexOnFace(e10[i], gf, pi);
fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 10);
}
for (unsigned int i = 0; i < e23.size(); i++){
SPoint2 pi; reparamMeshVertexOnFace(e23[i], gf, pi); fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 23);
}
for (unsigned int i = 0; i < e32.size(); i++){
SPoint2 pi; reparamMeshVertexOnFace(e32[i], gf, pi);
fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 32);
}
for (unsigned int i = 0; i < e02.size(); i++){
SPoint2 pi; reparamMeshVertexOnFace(e02[i], gf, pi);
fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 2);
}
for (unsigned int i = 0; i < e13.size(); i++){
SPoint2 pi; reparamMeshVertexOnFace(e13[i], gf, pi);
fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 13);
}
fprintf(f2,"};\n");
fclose(f2);
char name3[234];
sprintf(name3,"quadXYZ_%d.pos", gf->tag());
FILE *f3 = Fopen(name3,"w");
if(!f3){
Msg::Error("Could not open file '%s'", name2);
return;
}
fprintf(f3,"View \"%s\" {\n",name3);
for (unsigned int i = 0; i < quads.size(); i++){
quads[i]->writePOS(f3,true,false,false,false,false,false);
}
fprintf(f3,"};\n");
fclose(f3);
}
static void createQuadsFromUV(GFace *gf, fullMatrix<SPoint2> &uv,
std::vector<std::vector<MVertex*> > &tab,
std::vector<MQuadrangle*> &newq, std::vector<MVertex*> &newv)
{
newq.clear();
newv.clear();
int M = uv.size1();
int N = uv.size2();
for (int i = 1; i < M-1; i++){
for (int j = 0; j < N-1; j++){
GPoint gp = gf->point(uv(i,j));
if (!gp.succeeded()){
printf("** QUADS FROM UV new vertex not created p=%g %g \n", uv(i,j).x(), uv(i,j).y()); //exit(1);
}
tab[i][j] = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
}
}
for (int i = 1; i < M-1; i++) tab[i][N-1] = tab[i][0];
//create vertices
for (int i = 1; i < M-1; i++)
for (int j = 0; j < N-1; j++)
newv.push_back(tab[i][j]);
// create quads
for (int i=0;i<M-1;i++){
for (int j=0;j<N-1;j++){
MQuadrangle *qnew = new MQuadrangle (tab[i][j],tab[i][j+1],tab[i+1][j+1],tab[i+1][j]);
newq.push_back(qnew);
}
}
return;
}
static std::vector<MVertex*> saturateEdgeRegular (GFace *gf, SPoint2 p1, SPoint2 p2,
int M, std::vector<SPoint2> &pe)
{
std::vector<MVertex*> pts;
for (int i=1;i<M;i++){
double s = ((double)i/((double)(M)));
SPoint2 p = p1 + (p2-p1)*s;
pe.push_back(p);
GPoint gp = gf->point(p);
if (!gp.succeeded()) printf("** SATURATE EDGE new vertex not created p=%g %g \n", p.x(), p.y());
MVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
if (!v){ pts.clear(); pts.resize(0); return pts;}
pts.push_back(v);
}
return pts;
}
static std::vector<MVertex*> saturateEdgeHarmonic (GFace *gf, SPoint2 p1, SPoint2 p2,
double H, double L,
int M, std::vector<SPoint2> &pe)
{
std::vector<MVertex*> pts;
for (int i=1;i<M;i++){
double y = ((double)(i))*H/M;
double Yy = cosh(M_PI*y/L) - tanh(M_PI*H/L)*sinh(M_PI*y/L);
double YyH = cosh(M_PI*H/L) - tanh(M_PI*H/L)*sinh(M_PI*H/L);
double s = (1 - Yy)/(1.-YyH);
SPoint2 p = p1 + (p2-p1)*s;
pe.push_back(p);
GPoint gp = gf->point(p);
if (!gp.succeeded()) printf("** SATURATE EDGE new vertex not created p=%g %g \n", p.x(), p.y());
MVertex *v = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
if (!v){ pts.clear(); pts.resize(0); return pts;}
pts.push_back(v);
}
return pts;
}
static void transfiniteSmoother(GFace* gf,
fullMatrix<SPoint2> &uv,
std::vector<std::vector<MVertex*> > &tab,
std::vector<MQuadrangle*> &newq,
std::vector<MVertex*> &newv,
bool isPeriodic=false)
{
int M = uv.size1();
int N = uv.size2();
int jStart = isPeriodic ? 0 : 1;
int numSmooth = 150;
fullMatrix<SPoint2> uvold = uv;
for(int k = 0; k < numSmooth; k++){
double norm = 0.0;
for (int i=1; i<M-1; i++){
for (int j = jStart; j<N-1; j++){
int jm1 = (j==0) ? N-2: j-1;
int jp1 = (isPeriodic ) ? (j+1)%(N-1) : j+1;
double alpha = 0.25 * (SQU(uv(i,jp1).x() - uv(i,jm1).x()) +
SQU(uv(i,jp1).y() - uv(i,jm1).y())) ;
double gamma = 0.25 * (SQU(uv(i+1,j).x() - uv(i-1,j).x()) +
SQU(uv(i+1,j).y() - uv(i-1,j).y()));
double beta = 0.0625 * ((uv(i+1,j).x() - uv(i-1,j).x()) * (uv(i,jp1).x() - uv(i,jm1).x()) +
(uv(i+1,j).y() - uv(i-1,j).y()) * (uv(i,jp1).y() - uv(i,jm1).y()));
double uk = 0.5 *
(alpha * (uv(i+1,j).x() + uv(i-1,j).x()) +
gamma * (uv(i,jp1).x() + uv(i,jm1).x()) -
2. * beta * (uv(i+1,jp1).x() - uv(i-1,jp1).x() -
uv(i+1,jm1).x() + uv(i-1,jm1).x())) / (alpha + gamma);
double vk = 0.5 *
(alpha * (uv(i+1,j).y() + uv(i-1,j).y()) +
gamma * (uv(i,jp1).y()+ uv(i,jm1).y()) -
2. * beta * (uv(i+1,jp1).y() - uv(i-1,jp1).y() -
uv(i+1,jm1).y() + uv(i-1,jm1).y())) / (alpha + gamma);
uvold(i,j) = uv(i,j);
norm += sqrt((uv(i,j).x()-uk)*(uv(i,j).x()-uk)+(uv(i,j).y()-vk)*(uv(i,j).y()-vk));
uv(i,j) = SPoint2(uk,vk);
}
}
if (isPeriodic){
for (int i = 1; i<M-1; i++) {
uv(i,N-1) = uv(i,0);
uvold(i,N-1) = uvold(i,0);
}
}
if(k%20==0){
printf("Transfinite smoother iter (%d) = %g \n", k, norm);
createQuadsFromUV(gf, uv, tab, newq, newv);
printQuads(gf, uv, newq, k+1);
}
}
//final print
createQuadsFromUV(gf, uv, tab, newq, newv);
printQuads(gf, uv, newq, numSmooth);
}
//elliptic surface grid generator
//see eq. 9.26 page 9-24 in handbook of grid generation
static void ellipticSmoother( GFace* gf,
fullMatrix<SPoint2> &uv,
std::vector<std::vector<MVertex*> > &tab,
std::vector<MQuadrangle*> &newq,
std::vector<MVertex*> &newv,
bool isPeriodic=false)
{
printQuads(gf, uv, newq, 0);
int nbSmooth = 70;
int M = uv.size1();
int N = uv.size2();
int jStart = isPeriodic ? 0 : 1;
int jEnd = N-1;
fullMatrix<SPoint2> uvold = uv;
for (int k = 0; k< nbSmooth; k++){
double norm = 0.0;
for (int i=1; i<M-1; i++){
for (int j = jStart; j<jEnd; j++){
int jm1 = (j==0) ? N-2: j-1;
int jp1 = (isPeriodic ) ? (j+1)%(N-1) : j+1;
Pair<SVector3, SVector3> der = gf->firstDer(uv(i,j));
SVector3 xu = der.first();
SVector3 xv = der.second();
// double g11 = dot(xu,xu);
// double g12 = dot(xu,xv);
// double g22 = dot(xv,xv);
SVector3 xuu = SVector3();
SVector3 xvv = SVector3();
SVector3 xuv = SVector3(); gf->secondDer(uv(i,j), &xuu, &xvv, &xuv);
double g11_bar = dot(xu,xu);
double g12_bar = dot(xu,xv);
double g22_bar = dot(xv,xv);
double dxsi = 1.;
double deta = 1.;
double u_xsi = (uv(i,jp1).x()-uv(i,j).x())/dxsi;
double u_eta = (uv(i+1,j).x()-uv(i,j).x())/deta;
double v_xsi = (uv(i,jp1).y()-uv(i,j).y())/dxsi;
double v_eta = (uv(i+1,j).y()-uv(i,j).y())/deta;
double g11 = g11_bar*u_xsi*u_xsi+2.*g12_bar*u_xsi*v_xsi+g22_bar*v_xsi*v_xsi;
double g12 = g11_bar*u_xsi*u_eta+g12_bar*(u_xsi*v_eta+u_eta*v_xsi)+g22_bar*v_xsi*v_eta;
double g22 = g11_bar*u_eta*u_eta+2.*g12_bar*u_eta*v_eta+g22_bar*v_eta*v_eta;
double jac = u_xsi*v_eta-u_eta*v_xsi; //sqrt(g11*g22-g12*g12);
double jac_bar = sqrt(g11_bar*g22_bar-g12_bar*g12_bar);
double g11_bar_u = 2.*dot(xu,xuu);
double g11_bar_v = 2.*dot(xu,xuv);
double g22_bar_u = 2.*dot(xv,xuv);
double g22_bar_v = 2.*dot(xv,xvv);
double g12_bar_u = dot(xu,xuv)+dot(xv,xuu);
double g12_bar_v = dot(xu,xvv)+dot(xv,xuv);
double jac_bar_u = 1./(2.*jac_bar)*(g11_bar*g22_bar_u+g22_bar*g11_bar_u-2.*g12_bar*g12_bar_u);
double jac_bar_v = 1./(2.*jac_bar)*(g11_bar*g22_bar_v+g22_bar*g11_bar_v-2.*g12_bar*g12_bar_v);
double lapl_u = 1./jac_bar*(jac_bar*(g22_bar_u-g12_bar_v)-(g22_bar*jac_bar_u-g12_bar*jac_bar_v));
double lapl_v = 1./jac_bar*(jac_bar*(g11_bar_v-g12_bar_u)-(g11_bar*jac_bar_v-g12_bar*jac_bar_u));
double over = 1./(4.*(g11+g22));
double uk = over*(2.*g22*(uv(i+1,j).x()+uv(i-1,j).x())+
2.*g11*(uv(i,jp1).x()+uv(i,jm1).x())-
2*g12*(uv(i+1,jp1).x()-uv(i-1,jp1).x()-uv(i+1,jm1).x()+uv(i-1,jm1).x())-
2.*jac*jac*lapl_u);
double vk = over*(2.*g22*(uv(i+1,j).y()+uv(i-1,j).y())+
2.*g11*(uv(i,jp1).y()+uv(i,jm1).y())-
2.*g12*(uv(i+1,jp1).y()-uv(i-1,jp1).y()-uv(i+1,jm1).y()+uv(i-1,jm1).y())-
2.*jac*jac*lapl_v);
uvold(i,j) = uv(i,j);
norm += sqrt((uv(i,j).x()-uk)*(uv(i,j).x()-uk)+(uv(i,j).y()-vk)*(uv(i,j).y()-vk));
uv(i,j) = SPoint2(uk,vk);
}
}
if (isPeriodic){
for (int i = 1; i<M-1; i++) {
uv(i,N-1) = uv(i,0);
uvold(i,N-1) = uvold(i,0);
}
}
//under-relaxation
// double omega = 0.8;
// for (int i=0; i<M; i++){
// for (int j = 0; j<N; j++){
// uv(i,j) = SPoint2(omega*uv(i,j).x() + (1.-omega)*uvold(i,j).x(),
// omega*uv(i,j).y() + (1.-omega)*uvold(i,j).y());
// }
// }
if(k%10==0){
printf("Elliptic smoother iter (%d) = %g \n", k, norm);
createQuadsFromUV(gf, uv, tab, newq, newv);
printQuads(gf, uv, newq, k+1);
}
}
createQuadsFromUV(gf, uv, tab, newq, newv);
printQuads(gf, uv, newq, nbSmooth);
}
//create initial grid points MxN using transfinite interpolation
/* c4 N c3
+--------------------------+
| | | | |
| | | | |
+--------------------------+ M
| | | | |
| | | | |
+--------------------------+
c1 c2
N
*/
static void createRegularGrid (GFace *gf,
MVertex *v1, SPoint2 &c1,
std::vector<MVertex*> &e12, std::vector<SPoint2> &pe12, int sign12,
MVertex *v2, SPoint2 &c2,
std::vector<MVertex*> &e23, std::vector<SPoint2> &pe23, int sign23,
MVertex *v3, SPoint2 &c3,
std::vector<MVertex*> &e34,std::vector<SPoint2> &pe34, int sign34,
MVertex *v4, SPoint2 &c4,
std::vector<MVertex*> &e41, std::vector<SPoint2> &pe41,int sign41,
std::vector<MQuadrangle*> &q,
std::vector<MVertex*> &newv,
fullMatrix<SPoint2> &uv,
std::vector<std::vector<MVertex*> > &tab)
{
int M = e23.size();
int N = e12.size();
uv.resize(M+2,N+2);
//std::vector<std::vector<MVertex*> > tab(M+2);
tab.resize(M+2);
for(int i = 0; i <= M+1; i++) tab[i].resize(N+2);
tab[0][0] = v1; uv(0,0) = c1;
tab[0][N+1] = v2; uv(0,N+1) = c2;
tab[M+1][N+1] = v3; uv(M+1,N+1) = c3;
tab[M+1][0] = v4; uv(M+1,0) = c4;
for (int i=0;i<N;i++){
tab[0][i+1] = sign12 > 0 ? e12 [i] : e12 [N-i-1];
tab[M+1][i+1] = sign34 < 0 ? e34 [i] : e34 [N-i-1];
uv(0,i+1) = sign12 > 0 ? pe12 [i] : pe12 [N-i-1];
uv(M+1,i+1) = sign34 < 0 ? pe34 [i] : pe34 [N-i-1];
}
for (int i=0;i<M;i++){
tab[i+1][N+1] = sign23 > 0 ? e23 [i] : e23 [M-i-1];
tab[i+1][0] = sign41 < 0 ? e41 [i] : e41 [M-i-1];
uv(i+1,N+1) = sign23 > 0 ? pe23 [i] : pe23 [M-i-1];
uv(i+1,0) = sign41 < 0 ? pe41 [i] : pe41 [M-i-1];
}
for (int i=0;i<N;i++){
for (int j=0;j<M;j++){
double u = (double) (i+1)/ ((double)(N+1));
double v = (double) (j+1)/ ((double)(M+1));
SPoint2 p12 = (sign12 >0) ? pe12[i] : pe12 [N-1-i];
SPoint2 p23 = (sign23 >0) ? pe23[j] : pe23 [M-1-j];
SPoint2 p34 = (sign34 <0) ? pe34[i] : pe34 [N-1-i];
SPoint2 p41 = (sign41 <0) ? pe41[j] : pe41 [M-1-j];
double Up = TRAN_QUAD(p12.x(), p23.x(), p34.x(), p41.x(),
c1.x(),c2.x(),c3.x(),c4.x(),u,v);
double Vp = TRAN_QUAD(p12.y(), p23.y(), p34.y(), p41.y(),
c1.y(),c2.y(),c3.y(),c4.y(),u,v);
if ((p12.x() && p12.y() == -1.0) || (p23.x() && p23.y() == -1.0) ||
(p34.x() && p34.y() == -1.0) || (p41.x() && p41.y() == -1.0)) {
Msg::Error("Wrong param -1");
return;
}
SPoint2 pij(Up,Vp);
GPoint gp = gf->point(pij);
if (!gp.succeeded()) printf("** INITIAL GRID new vertex not created p=%g %g \n", Up, Vp);
MVertex *vnew = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
newv.push_back(vnew);
uv(j+1,i+1) = pij;
tab[j+1][i+1] = vnew;
}
}
// create quads
for (int i=0;i<=N;i++){
for (int j=0;j<=M;j++){
MQuadrangle *qnew = new MQuadrangle (tab[j][i],tab[j][i+1],tab[j+1][i+1],tab[j+1][i]);
q.push_back(qnew);
}
}
}
//create initial grid points MxN using transfinite interpolation
/*
c2 N c2
+--------------------------+
| | | | |
| | | | |
+--------------------------+ M
| | | | |
| | | | |
+--------------------------+
c0 c0
N
*/
static void createRegularGridPeriodic (GFace *gf,int sign2,
MVertex *v0, SPoint2 &c0,
MVertex *v2, SPoint2 &c2,
std::vector<MVertex*> &e02, std::vector<SPoint2> &pe02,
std::vector<MVertex*> &e00, std::vector<SPoint2> &pe00,
std::vector<MVertex*> &e22, std::vector<SPoint2> &pe22,
std::vector<MQuadrangle*> &q,
std::vector<MVertex*> &newv,
fullMatrix<SPoint2> &uv,
std::vector<std::vector<MVertex*> > &tab)
{
int M = e02.size();
int N = e00.size();
uv.resize(M+2,N+2);
char name3[234];
sprintf(name3,"quadParam_%d.pos", gf->tag());
FILE *f3 = Fopen(name3,"w");
if(f3) fprintf(f3,"View \"%s\" {\n",name3);
tab.resize(M+2);
for(int i = 0; i < M+2; i++) tab[i].resize(N+2);
//periodic boundary mesh vertices
tab[0][0] = v0; uv(0,0) = c0;
tab[0][N+1] = v0; uv(0,N+1) = c0;
tab[M+1][N+1] = v2; uv(M+1,N+1) = c2;
tab[M+1][0] = v2; uv(M+1,0) = c2;
for (int i=0;i<N;i++){
tab[0][i+1] = e00 [i]; uv(0,i+1) = pe00 [i] ;
tab[M+1][i+1] = (sign2 > 0) ? e22 [i] : e22 [N-i-1] ;
uv(M+1,i+1) = (sign2 > 0) ? pe22 [i] : pe22 [N-i-1] ; }
for (int i=0;i<M;i++){
tab[i+1][N+1] = e02 [i]; uv(i+1,N+1) = pe02 [i] ;
tab[i+1][0] = e02 [i]; uv(i+1,0) = pe02 [i];
}
//interior mesh_vertices
for (int i=0;i<N;i++){
SPoint2 p0 = pe00[i] ;
SPoint2 p2 = (sign2>0) ? pe22[i] : pe22 [N-i-1] ;
std::vector<SPoint2> pei;
std::vector<MVertex*> ei = saturateEdgeRegular(gf,p0,p2,M+1,pei);//M+1
for (int j=0;j<M;j++){
SPoint2 pij = pei[j];
GPoint gp = gf->point(pij);
if (!gp.succeeded()) printf("** INITIAL GRID new vertex not created p=%g %g \n", pij.x(), pij.y());
tab[j+1][i+1] = new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
uv(j+1,i+1) = pij;
}
}
//create vertices
for (int i=0;i<N+1;i++)
for (int j=1;j<M+1;j++)
newv.push_back(tab[j][i]);
// create quads
for (int i=0;i<=N;i++){
for (int j=0;j<=M;j++){
MQuadrangle *qnew = new MQuadrangle (tab[j][i],tab[j][i+1],tab[j+1][i+1],tab[j+1][i]);
q.push_back(qnew);
}
}
//print
if(f3){
for (int i=0;i<N+2;i++)
for (int j=0;j<M+2;j++)
fprintf(f3,"SP(%g,%g,%g) {%d};\n", uv(j,i).x(), uv(j,i).y(), 0.0, j);
fprintf(f3,"};\n");
fclose(f3);
}
}
static void updateFaceQuads(GFace *gf, std::vector<MQuadrangle*> &quads, std::vector<MVertex*> &newv)
{
for (unsigned int i = 0; i < quads.size(); i++){
gf->quadrangles.push_back(quads[i]);
}
for(unsigned int i = 0; i < newv.size(); i++){
gf->mesh_vertices.push_back(newv[i]);
}
}
static bool computeRingVertices(GFace *gf, Centerline *center,
std::vector<MVertex*> &vert1, std::vector<MVertex*> &vert2,
int &N, int &M, int &close_ind, int &sign2,
double &arc, double &length)
{
#if defined(HAVE_ANN)
std::list<GEdge*> bedges = gf->edges();
std::list<GEdge*>::iterator itb = bedges.begin();
std::list<GEdge*>::iterator ite = bedges.end(); ite--;
GEdge *ge1 = (*itb)->getCompound();
GEdge *ge2 = (*ite)->getCompound();
N = ge1->mesh_vertices.size() + 1;
int N2 = ge2->mesh_vertices.size() + 1;
if (N != N2 || N%2 != 0){
Msg::Error("You should an even number nbPoints in centerline =%d \n", N); exit(1);
return false;
}
vert1.push_back(ge1->getBeginVertex()->mesh_vertices[0]);
vert2.push_back(ge2->getBeginVertex()->mesh_vertices[0]);
for(int i = 0; i < N-1; i++) {
vert1.push_back(ge1->mesh_vertices[i]);
vert2.push_back(ge2->mesh_vertices[i]);
}
SPoint2 pv1; reparamMeshVertexOnFace(vert1[0], gf, pv1);
SPoint2 pv2; reparamMeshVertexOnFace(vert2[0], gf, pv2);
SPoint2 pv1b; reparamMeshVertexOnFace(vert1[1], gf, pv1b);
SPoint2 pv2b; reparamMeshVertexOnFace(vert2[1], gf, pv2b);
SPoint2 pv1c; reparamMeshVertexOnFace(vert1[2], gf, pv1c);
SPoint2 pv2c; reparamMeshVertexOnFace(vert2[2], gf, pv2c);
SVector3 vec1(pv1b.x()-pv1.x(),pv1b.y()-pv1.y() , 0.0);
SVector3 vec1b(pv1c.x()-pv1b.x(),pv1c.y()-pv1b.y() ,0.0);
SVector3 vec2(pv2b.x()-pv2.x(),pv2b.y()-pv2.y() ,0.0);
SVector3 vec2b(pv2c.x()-pv2b.x(),pv2c.y()-pv2b.y() , 0.0);
sign2 = (dot(crossprod(vec1,vec1b),crossprod(vec2,vec2b)) < 0) ? -1 : +1 ;
double n1 = sqrt(pv1.x()*pv1.x()+pv1.y()*pv1.y());
double n2 = sqrt(pv2.x()*pv2.x()+pv2.y()*pv2.y());
std::vector<MVertex*> vert_temp;
if (n2 > n1) {
vert_temp = vert1;
vert1.clear();
vert1 = vert2;
vert2.clear();
vert2 = vert_temp;
}
double rad = center->operator()(vert1[0]->x(), vert1[0]->y(), vert1[0]->z());
ANNpointArray nodes = annAllocPts(N, 3);
ANNidxArray index = new ANNidx[1];
ANNdistArray dist = new ANNdist[1];
for (int ind = 0; ind < N; ind++){
SPoint2 pp2; reparamMeshVertexOnFace(vert2[ind], gf, pp2);
nodes[ind][0] = pp2.x();
nodes[ind][1] = pp2.y();
nodes[ind][2] = 0.0;
}
ANNkd_tree *kdtree = new ANNkd_tree(nodes, N, 3);
SPoint2 pp1; reparamMeshVertexOnFace(vert1[0], gf, pp1);
double xyz[3] = {pp1.x(), pp1.y(),0.0};
kdtree->annkSearch(xyz, 1, index, dist);
close_ind = index[0];
SPoint2 p1; reparamMeshVertexOnFace(vert1[0], gf, p1);
SPoint2 p2; reparamMeshVertexOnFace(vert2[close_ind], gf, p2);
GPoint gp_i(vert1[0]->x(), vert1[0]->y(), vert1[0]->z());
SPoint2 p1b; reparamMeshVertexOnFace(vert1[N/2], gf, p1b);
SPoint2 p2b; reparamMeshVertexOnFace(vert2[(close_ind+N/2)%N], gf, p2b);
GPoint gp_ib(vert1[N/2]->x(), vert1[N/2]->y(), vert1[N/2]->z());
length = 0.0;
double lengthb = 0.0;
int nb = 100;
for (int i = 0; i< nb; i++){
SPoint2 pii(p1.x() + (double)(i+1)/nb*(p2.x()-p1.x()), p1.y() + (double)(i+1)/nb*(p2.y()-p1.y()));
GPoint gp_ii = gf->point(pii);
SPoint2 piib(p1b.x() + (double)(i+1)/nb*(p2b.x()-p1b.x()), p1b.y() + (double)(i+1)/nb*(p2b.y()-p1b.y()));
GPoint gp_iib = gf->point(piib);
length += gp_i.distance(gp_ii);
lengthb += gp_ib.distance(gp_iib);
gp_i = gp_ii;
gp_ib = gp_iib;
}
arc = 2*M_PI*rad;
double lc = arc/N; M = (int)(0.5*(length+lengthb)/lc) ;
delete kdtree;
delete[]index;
delete[]dist;
annDeallocPts(nodes);
return true;
#else
return false;
#endif
}
//vert1 is the outer circle
//vert2 is the inner circle
// - - - -
// - -
// - - -
// v0- v2- - -
// - - -
// - -
// - - - -
bool createRegularTwoCircleGridPeriodic (Centerline *center, GFace *gf)
{
#if defined(HAVE_ANN)
std::vector<MVertex*> vert1, vert2;
int N, M, close_ind, sign2;
double arc, length;
bool success = computeRingVertices(gf, center, vert1, vert2, N, M, close_ind, sign2, arc,length);
if(!success) return false;
MVertex *v0 = vert1[0]; SPoint2 p0; reparamMeshVertexOnFace(v0, gf, p0);
MVertex *v2 = vert2[close_ind]; SPoint2 p2; reparamMeshVertexOnFace(v2, gf, p2);
//printf("grid N = %d M = %d\n", N , M);
std::vector<MVertex*> e00,e22;//edges without first and last vertex
for (int i=1;i<N;i++) e00.push_back(vert1[i]);
for (int i=1;i<N;i++) e22.push_back(vert2[(close_ind+i)%N]);
std::vector<SPoint2> pe00, pe22;
for (unsigned int i = 0; i < e00.size(); i++){
SPoint2 p00; reparamMeshVertexOnFace(e00[i], gf, p00);
SPoint2 p22; reparamMeshVertexOnFace(e22[i], gf, p22);
pe00.push_back(p00);
pe22.push_back(p22);
}
std::vector<SPoint2> pe02;
std::vector<MVertex*> e02 = saturateEdgeRegular(gf,p0,p2,M+1,pe02);
std::vector<MQuadrangle*> quads;
std::vector<MVertex*> newv;
fullMatrix<SPoint2> uv;
std::vector<std::vector<MVertex*> > tab;
createRegularGridPeriodic (gf,sign2,
v0, p0,
v2, p2,
e02, pe02,
e00, pe00,
e22, pe22,
quads, newv, uv, tab);
printQuads(gf, uv, quads, 0);
transfiniteSmoother(gf, uv, tab, quads, newv, true);
//ellipticSmoother(gf, uv, tab, quads, newv, true);
updateFaceQuads(gf, quads, newv);
//exit(1);
//printParamGrid(gf, vert1, vert2, e00,e22,e02,e02,e02,e02, quads);
return true;
#else
return false;
#endif
}
//vert1 is the outer circle
//vert2 is the inner circle
// - - - -
// - -
// - - -
// v0- v2- -v3 -v1
// - - -
// - -
// - - - -
bool createRegularTwoCircleGrid (Centerline *center, GFace *gf)
{
#if defined(HAVE_ANN)
std::vector<MVertex*> vert1, vert2;
int N, M, close_ind, sign2;
double arc, length;
bool success = computeRingVertices(gf, center, vert1, vert2, N, M, close_ind, sign2, arc,length);
if(!success) return false;
MVertex *v0 = vert1[0]; SPoint2 p0; reparamMeshVertexOnFace(v0, gf, p0);
MVertex *v1 = vert1[N/2]; SPoint2 p1; reparamMeshVertexOnFace(v1, gf, p1);
MVertex *v2 = vert2[close_ind]; SPoint2 p2; reparamMeshVertexOnFace(v2, gf, p2);
MVertex *v3 = vert2[(close_ind+N/2)%N]; SPoint2 p3; reparamMeshVertexOnFace(v3, gf, p3);
printf("grid N = %d M = %d\n", N , M);
std::vector<MVertex*> e01,e10,e23,e32;//edges without first and last vertex
for (int i=1;i<N/2;i++) e01.push_back(vert1[i]);
for (int i=N/2+1;i<N;i++) e10.push_back(vert1[i]);
for (int i=1;i<N/2;i++) e23.push_back(vert2[(close_ind+i)%N]);
for (int i=N/2+1;i<N;i++) e32.push_back(vert2[(close_ind+i)%N]);
std::vector<SPoint2> pe01, pe10, pe23, pe32;
for (unsigned int i = 0; i < e01.size(); i++){
SPoint2 p01; reparamMeshVertexOnFace(e01[i], gf, p01);
SPoint2 p10; reparamMeshVertexOnFace(e10[i], gf, p10);
SPoint2 p23; reparamMeshVertexOnFace(e23[i], gf, p23);
SPoint2 p32; reparamMeshVertexOnFace(e32[i], gf, p32);
pe01.push_back(p01);
pe10.push_back(p10);
pe23.push_back(p23);
pe32.push_back(p32);
}
std::vector<SPoint2> pe02, pe13;
std::vector<MVertex*> e02 = saturateEdgeHarmonic (gf,p0,p2,length, arc, M+1, pe02);
std::vector<MVertex*> e13 = saturateEdgeHarmonic (gf,p1,p3,length, arc, M+1, pe13);
std::vector<MVertex*> e_inner1 = e23;
std::vector<MVertex*> e_inner2 = e32;
std::vector<SPoint2> pe_inner1 = pe23;
std::vector<SPoint2> pe_inner2 = pe32;
if (sign2 <0){
e_inner1 = e32;
e_inner2 = e23;
pe_inner1 = pe32;
pe_inner2 = pe23;
}
std::vector<MQuadrangle*> quads, quadS1, quadS2;
std::vector<MVertex*> newv, newv1, newv2;
fullMatrix<SPoint2> uv1, uv2; std::vector<std::vector<MVertex*> > tab1, tab2;
createRegularGrid (gf,
v0, p0,
e01, pe01, +1,
v1,p1,
e13, pe13, +1,
v3,p3,
e_inner1, pe_inner1, -sign2,
v2,p2,
e02, pe02, -1,
quads, newv, uv1, tab1);
createRegularGrid (gf,
v0,p0,
e02, pe02, +1,
v2,p2,
e_inner2, pe_inner2, -sign2,
v3,p3,
e13, pe13, -1,
v1,p1,
e10, pe10, +1,
quads, newv, uv2, tab2);
printQuads(gf, uv1, quads, 0);
ellipticSmoother(gf, uv1, tab1, quadS1, newv1);
ellipticSmoother(gf, uv2, tab2, quadS2, newv2);
quads.clear();
for (unsigned int i= 0; i< quadS1.size(); i++) quads.push_back(quadS1[i]);
for (unsigned int i= 0; i< quadS2.size(); i++) quads.push_back(quadS2[i]);
//updateFaceQuads(gf, quads, newv);
printParamGrid(gf, vert1, vert2, e01,e10,e23,e32,e02,e13, quads);
return true;
#else
return false;
#endif
}