Skip to content
Snippets Groups Projects
Commit 67723cda authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

bug fix in the 2D mesher

parent 7fca694b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/SBoundingBox3d.h
+ 7
− 8
View file @ 67723cda
......@@ -8,6 +8,7 @@
#include <float.h>
#include "SPoint3.h"
#include "SVector3.h"
#if defined(WIN32)
#undef min
......@@ -86,14 +87,12 @@ class SBoundingBox3d {
SPoint3 center() const { return (MinPt + MaxPt) * .5; }
void makeCube()
{
SPoint3 len = MaxPt - MinPt;
double scales[3];
double max = -1.0;
for(int i = 0; i < 3; i++)
max = len[i] > max ? len[i] : max;
for(int j = 0; j < 3; j++)
scales[j] = max/len[j];
scale(scales[0],scales[1],scales[2]);
SVector3 len = MaxPt - MinPt;
SPoint3 cc = center();
MaxPt = cc + SPoint3(1,1,1);
MinPt = cc + SPoint3(-1,-1,-1);
double sc = len.norm() * 0.5;
scale (sc,sc,sc);
}
private:
SPoint3 MinPt, MaxPt;
......
Mesh/meshGFace.cpp
+ 37
− 0
View file @ 67723cda
......@@ -705,6 +705,17 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
gf->meshStatistics.status = GFace::DONE;
return true;
}
if(all_vertices.size() == 3){
MVertex *vv[3];
int i = 0;
for(std::set<MVertex*>::iterator it = all_vertices.begin();
it != all_vertices.end(); it++){
vv[i++] = *it;
}
gf->triangles.push_back(new MTriangle(vv[0], vv[1], vv[2]));
gf->meshStatistics.status = GFace::DONE;
return true;
}
// Buid a BDS_Mesh structure that is convenient for doing the actual
// meshing procedure
......@@ -735,6 +746,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// here check if some boundary layer nodes should be added
bbox.makeCube();
// compute the bounding box in parametric space
SVector3 dd(bbox.max(), bbox.min());
double LC2D = norm(dd);
......@@ -1456,6 +1469,26 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
if(nbPointsTotal < 3){
Msg::Warning("Mesh Generation of Model Face %d Skipped: "
"Only %d Mesh Vertices on The Contours",
gf->tag(), nbPointsTotal);
gf->meshStatistics.status = GFace::DONE;
return true;
}
if(nbPointsTotal == 3){
MVertex *vv[3];
int i = 0;
for(std::map<BDS_Point*, MVertex*>::iterator it = recoverMap.begin();
it != recoverMap.end(); it++){
vv[i++] = it->second;
}
gf->triangles.push_back(new MTriangle(vv[0], vv[1], vv[2]));
gf->meshStatistics.status = GFace::DONE;
return true;
}
// Use a divide & conquer type algorithm to create a triangulation.
// We add to the triangulation a box with 4 points that encloses the
// domain.
......@@ -1481,6 +1514,10 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// Increase the size of the bounding box, add 4 points that enclose
// the domain, use negative number to distinguish those fake
// vertices
// FIX A BUG HERE IF THE SIZE OF THE BOX IS ZERO
bbox.makeCube();
bbox *= 3.5;
MVertex *bb[4];
bb[0] = new MVertex(bbox.min().x(), bbox.min().y(), 0, 0, -1);
......
Mesh/meshGFaceDelaunayInsertion.cpp
+ 72
− 97
View file @ 67723cda
......@@ -1054,6 +1054,7 @@ double optimalPointFrontal (GFace *gf,
// printf("%12.5E %12.5E\n",d,RATIO);
// const double L = d ;
// avoid to go toooooo far
const double L = d > q ? q : d;
newPoint[0] = midpoint[0] + L * dir[0]/RATIO;
......@@ -1215,6 +1216,7 @@ void bowyerWatsonFrontal(GFace *gf)
#endif
}
void optimalPointFrontalQuad (GFace *gf,
MTri3* worst,
int active_edge,
......@@ -1298,15 +1300,20 @@ void optimalPointFrontalQuad (GFace *gf,
}
}
void optimalPointFrontalQuadAniso (GFace *gf,
MTri3* worst,
int active_edge,
std::vector<double> &Us,
std::vector<double> &Vs,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
double newPoint[2],
double metric[3]){
void optimalPointFrontalQuadB (GFace *gf,
MTri3* worst,
int active_edge,
std::vector<double> &Us,
std::vector<double> &Vs,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
double newPoint[2],
double metric[3]){
// optimalPointFrontalQuad (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
MTriangle *base = worst->tri();
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
......@@ -1328,100 +1335,67 @@ void optimalPointFrontalQuadAniso (GFace *gf,
// rotate the points with respect to the angle
double XP1 = 0.5*(Q[0] - P[0]);
double YP1 = 0.5*(Q[1] - P[1]);
double DX = 0.5*(Q[0] + P[0]);
double DY = 0.5*(Q[0] + P[0]);
double xp = XP1 * cos(quadAngle) + YP1 * sin(quadAngle);
double yp = -XP1 * sin(quadAngle) + YP1 * cos(quadAngle);
double X4 = O[0] -DX;
double Y4 = O[1] -DY;
double x4 = X4 * cos(quadAngle) + Y4 * sin(quadAngle);
double y4 = -X4 * sin(quadAngle) + Y4 * cos(quadAngle);
// ensure xp > yp
bool exchange = false;
if (fabs(xp) < fabs(yp)){
double temp = xp;
xp = yp;
yp = temp;
temp = x4;
x4 = y4;
y4 = temp;
exchange = true;
double t = 0.0;
if (fabs(xp)>=fabs(yp)){
t = -1. + 2*fabs(xp*yp)/(xp*xp + yp*yp);
}
buildMetric(gf, midpoint, metric);
double RATIO = 1./pow(metric[0]*metric[2]-metric[1]*metric[1],0.25);
const double rhoM1 = 0.5 * RATIO *
else {
t = 1. - 2*fabs(xp*yp)/(xp*xp + yp*yp);
}
const double usq = 1./sqrt(2.0);
double factor = usq + fabs(t) * (1.-usq);
// printf("t = %lf %lf\n",t,factor);
SVector3 P3(base->getVertex(ip1)->x(),base->getVertex(ip1)->y(),base->getVertex(ip1)->z());
SVector3 Q3(base->getVertex(ip2)->x(),base->getVertex(ip2)->y(),base->getVertex(ip2)->z());
SVector3 O3(base->getVertex(ip3)->x(),base->getVertex(ip3)->y(),base->getVertex(ip3)->z());
SVector3 PMID = P3*(1.-t)*.5 + Q3*(t+1.)*.5;
SVector3 N = crossprod (O3-P3,Q3-P3);
SVector3 T = Q3-P3;
N.normalize();
T.normalize();
SVector3 N2 = crossprod (T,N);
N2.normalize();
if (dot (N2,O3-PMID) < 0) N2 = N2*(-1.0);
const double rhoM1 = 0.5 *
(vSizes[base->getVertex(ip1)->getIndex()] +
vSizes[base->getVertex(ip2)->getIndex()] );
const double rhoM2 = 0.5 * RATIO *
vSizes[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
const double rhoM2 = 0.5 *
(vSizesBGM[base->getVertex(ip1)->getIndex()] +
vSizesBGM[base->getVertex(ip2)->getIndex()] );// * RATIO;
const double rhoM = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
double npx,npy;
const double L_edge = lengthInfniteNorm(P,Q,quadAngle);
const double XP[2]={fabs(xp),fabs(yp)};
if (xp*yp > 0){
double xe[2] = { rhoM - fabs(xp) + fabs(yp), rhoM };
double xc[2] = {0.5*(x4 - fabs(xp)), 0.5*(x4 + fabs(xp)) - fabs(yp)};
double xl[2] = {rhoM, rhoM+fabs(xp)-fabs(yp)};
const double R_Active = lengthInfniteNorm(xc,XP,0.0); // alerady rotated !
const double L_ec = lengthInfniteNorm(xe,xc,0.0); // alerady rotated !
const double L_el = lengthInfniteNorm(xe,xl,0.0); // alerady rotated !
double t = ( L_edge - rhoM)/(L_edge - R_Active);
t = std::max(1.,t);
t = std::min(-L_el/L_ec,t);
npx = xe[0] + t*(xc[0]-xe[0]);
npy = xe[1] + t*(xc[1]-xe[1]);
}
else {
double xe[2] = { rhoM - xp + yp, rhoM };
double xc[2] = {0.5*(x4 - xp), 0.5*(x4 + xp) - yp};
double xl[2] = {rhoM, rhoM+xp-yp};
const double R_Active = lengthInfniteNorm(xc,XP,0.0); // alerady rotated !
const double L_ec = lengthInfniteNorm(xe,xc,0.0); // alerady rotated !
const double L_el = lengthInfniteNorm(xe,xl,0.0); // alerady rotated !
const double XP[2]={xp,yp};
double t = ( L_edge - rhoM)/(L_edge - R_Active);
t = std::max(1.,t);
t = std::min(-L_el/L_ec,t);
npx = xe[0] + t*(xc[0]-xe[0]);
npy = xe[1] + t*(xc[1]-xe[1]);
vSizesBGM[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
const double a = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
double d = a*factor;
if (gf->geomType() == GEntity::CompoundSurface){
GFaceCompound *gfc = dynamic_cast<GFaceCompound*> (gf);
if (gfc){
GPoint gp = gfc->intersectionWithCircle(N2,N,PMID,d,newPoint);
if (gp.succeeded()){
newPoint[0] = gp.u();
newPoint[1] = gp.v();
return ;
}
}
}
double uvt[3] = {newPoint[0],newPoint[1],0.0};
curveFunctorCircle cc (N2,N,PMID,d);
surfaceFunctorGFace ss (gf);
/*
if (xp*yp > 0){
npx = - fabs(xp)*factor;
npy = fabs(xp)*(1.+factor) - fabs(yp);
if (intersectCurveSurface (cc,ss,uvt,d*1.e-8)){
// printf("%g %g vs %g %g\n", newPoint[0], newPoint[1],uvt[0],uvt[1]);
newPoint[0] = uvt[0];
newPoint[1] = uvt[1];
}
else {
npx = fabs(xp) * factor;
npy = (1.+factor)*fabs(xp) - fabs(yp);
}
*/
if (exchange){
double temp = npx;
npx = npy;
npy = temp;
}
newPoint[0] = midpoint[0] + cos(quadAngle) * npx - sin(quadAngle) * npy;
newPoint[1] = midpoint[1] + sin(quadAngle) * npx + cos(quadAngle) * npy;
if ((midpoint[0] - newPoint[0])*(midpoint[0] - O[0]) +
(midpoint[1] - newPoint[1])*(midpoint[1] - O[1]) < 0){
newPoint[0] = midpoint[0] - cos(quadAngle) * npx + sin(quadAngle) * npy;
newPoint[1] = midpoint[1] - sin(quadAngle) * npx - cos(quadAngle) * npy;
newPoint[0] = -10000;
newPoint[1] = 100000;
Msg::Warning("--- Non optimal point found -----------");
}
buildMetric(gf, newPoint, metric);
}
......@@ -1527,11 +1501,12 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
if (!ActiveTris.size())break;
// if (gf->tag() == 1900){ char name[245];
// sprintf(name,"x_GFace_%d_Layer_%d.pos",gf->tag(),ITER);
// _printTris (name, AllTris, Us,Vs,true);
// }
/* if (1 || gf->tag() == 1900){
char name[245];
sprintf(name,"x_GFace_%d_Layer_%d.pos",gf->tag(),ITER);
_printTris (name, AllTris, Us,Vs,true);
}
*/
std::set<MTri3*,compareTri3Ptr>::iterator WORST_ITER = ActiveTris.begin();
MTri3 *worst = (*WORST_ITER);
......
benchmarks/2d/square.geo
+ 6
− 17
View file @ 67723cda
//Field[1] = Attractor;
//Field[1].EdgesList = {1};
//Field[1].NNodesByEdge = 10;
//Background Field = 1;
//Field[1] = MathEval;
//Field[1].F = "1.0"; //0.1*x+0.1";
//Background Field = 1;
Mesh.Algorithm=1;
Mesh.CharacteristicLengthFactor=1.5;
lc=0.1;
Point(1) = {0, 0, 0}; //,lc};
Point(2) = {0, 10, 0}; //,lc};
Point(3) = {10, 10, 0}; //,lc};
Point(4) = {10, 0, 0}; //,lc};
lc=1;
Point(1) = {0, 0, 0,lc*.1};
Point(2) = {0, 10, 0,lc};
Point(3) = {10, 10, 0,lc};
Point(4) = {10, 0, 0,lc};
Line(1) = {2, 3};
Line(2) = {3, 4};
Line(3) = {4, 1};
......@@ -36,3 +24,4 @@ Plane Surface(10) = {5};
Recombine Surface {10};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment