diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 0fd1b5c576fe62d22cecd5b76bcdc7bdddc188b8..cf34e971fd8e413432d511760a3c71698bbb986f 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.109 2007-12-03 15:17:39 remacle Exp $
+// $Id: CommandLine.cpp,v 1.110 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -241,6 +241,10 @@ void Get_Options(int argc, char *argv[])
CTX.mesh.optimize = 1;
i++;
}
+ else if(!strcmp(argv[i] + 1, "optimize_netgen")) {
+ CTX.mesh.optimizeNetgen = 1;
+ i++;
+ }
else if(!strcmp(argv[i] + 1, "nopopup")) {
CTX.nopopup = 1;
i++;
@@ -353,6 +357,17 @@ void Get_Options(int argc, char *argv[])
exit(1);
}
}
+ }
+ else if(!strcmp(argv[i] + 1, "clmin")) {
+ i++;
+ if(argv[i] != NULL) {
+ CTX.mesh.lc_min = atof(argv[i++]);
+ if(CTX.mesh.lc_factor <= 0.0) {
+ fprintf(stderr, ERROR_STR
+ "Minimum length size must be > 0\n");
+ exit(1);
+ }
+ }
else {
fprintf(stderr, ERROR_STR "Missing number\n");
exit(1);
diff --git a/Common/Context.h b/Common/Context.h
index 7933a761dc2ae866d22daab443ee3b8378ae53d8..d5a5390c3bcd09c884759f4959ba23fe510b5e78 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -181,7 +181,7 @@ public :
int optimize, optimizeNetgen, refine_steps;
int quality_type, label_type;
double quality_inf, quality_sup, radius_inf, radius_sup;
- double scaling_factor, lc_factor, rand_factor, lc_integration_precision;
+ double scaling_factor, lc_factor, rand_factor, lc_integration_precision,lc_min;
int dual, draw_skin_only;
int light, light_two_side, light_lines;
int format, nb_smoothing, algo2d, algo3d, algo_recombine;
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index d12f23bda8cd97fad3535d635f199eb6a1a708a5..3d206e35d97fdaf3c49ed84633f9d6bdfa2e6c7f 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -895,6 +895,8 @@ StringXNumber MeshOptions_Number[] = {
"Field format for Nastran BDF files (0=free, 1=small, 2=large)" },
{ F|O, "CharacteristicLengthFactor" , opt_mesh_lc_factor , 1.0 ,
"Factor applied to all characteristic lengths" },
+ { F|O, "CharacteristicLengthFactor" , opt_mesh_lc_min, 0.0 ,
+ "Minimum mesh size" },
{ F|O, "CharacteristicLengthFromCurvature" , opt_mesh_lc_from_curvature , 0. ,
"Compute characteritic lenghts automatically from curvatures" },
{ F|O, "ColorCarousel" , opt_mesh_color_carousel , 1. ,
diff --git a/Common/Options.cpp b/Common/Options.cpp
index 0c64aa4d900f0ccb2bcaed2ddc74bf6ae88897fe..6622d40d37ab22e88258aec3928a44cc8789e02f 100644
--- a/Common/Options.cpp
+++ b/Common/Options.cpp
@@ -1,4 +1,4 @@
-// $Id: Options.cpp,v 1.374 2008-01-12 18:40:14 geuzaine Exp $
+// $Id: Options.cpp,v 1.375 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -4294,6 +4294,17 @@ double opt_mesh_lc_factor(OPT_ARGS_NUM)
return CTX.mesh.lc_factor;
}
+double opt_mesh_lc_min(OPT_ARGS_NUM)
+{
+ if(action & GMSH_SET)
+ CTX.mesh.lc_min = val;
+#if defined(HAVE_FLTK)
+ if(WID && (action & GMSH_GUI))
+ WID->mesh_value[25]->value(CTX.mesh.lc_min);
+#endif
+ return CTX.mesh.lc_min;
+}
+
double opt_mesh_lc_from_curvature(OPT_ARGS_NUM)
{
if(action & GMSH_SET)
diff --git a/Common/Options.h b/Common/Options.h
index 0de824f8824e442cfcf6c6121e21a48aa7343fdb..a507131e2b11024eebe666fb4733c49e8b5bcac5 100644
--- a/Common/Options.h
+++ b/Common/Options.h
@@ -433,6 +433,7 @@ double opt_mesh_normals(OPT_ARGS_NUM);
double opt_mesh_tangents(OPT_ARGS_NUM);
double opt_mesh_explode(OPT_ARGS_NUM);
double opt_mesh_scaling_factor(OPT_ARGS_NUM);
+double opt_mesh_lc_min(OPT_ARGS_NUM);
double opt_mesh_lc_factor(OPT_ARGS_NUM);
double opt_mesh_lc_from_curvature(OPT_ARGS_NUM);
double opt_mesh_rand_factor(OPT_ARGS_NUM);
diff --git a/Fltk/Callbacks.cpp b/Fltk/Callbacks.cpp
index 32f1b67047071f1b98581509656b7f49dd2cb791..2528340f2e2565f4f8916a919c2338682d0e0c4d 100644
--- a/Fltk/Callbacks.cpp
+++ b/Fltk/Callbacks.cpp
@@ -1,4 +1,4 @@
-// $Id: Callbacks.cpp,v 1.558 2008-01-12 18:40:14 geuzaine Exp $
+// $Id: Callbacks.cpp,v 1.559 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -1130,6 +1130,7 @@ void mesh_options_ok_cb(CALLBACK_ARGS)
opt_mesh_nb_smoothing(0, GMSH_SET, WID->mesh_value[0]->value());
opt_mesh_lc_factor(0, GMSH_SET, WID->mesh_value[2]->value());
+ opt_mesh_lc_min(0, GMSH_SET, WID->mesh_value[25]->value());
opt_mesh_quality_inf(0, GMSH_SET, WID->mesh_value[4]->value());
opt_mesh_quality_sup(0, GMSH_SET, WID->mesh_value[5]->value());
opt_mesh_radius_inf(0, GMSH_SET, WID->mesh_value[6]->value());
diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 2a10928305f778702dc6458ca38a7f6c582eb7f8..af5ab739c1c3e2137948fb937f58cc13e591f219 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.648 2008-01-12 18:40:14 geuzaine Exp $
+// $Id: GUI.cpp,v 1.649 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -2432,7 +2432,11 @@ void GUI::create_option_window()
mesh_value[2]->align(FL_ALIGN_RIGHT);
mesh_value[2]->callback(mesh_options_ok_cb);
- mesh_value[3] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 6 * BH, IW, BH, "Element order");
+ mesh_value[25] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 6 * BH, IW, BH, "Minimum mesh size");
+ mesh_value[25]->align(FL_ALIGN_RIGHT);
+ mesh_value[25]->callback(mesh_options_ok_cb);
+
+ mesh_value[3] = new Fl_Value_Input(L + 2 * WB, 2 * WB + 7 * BH, IW, BH, "Element order");
mesh_value[3]->minimum(1);
// FIXME: this makes it possible to set > 2 by hand, but not by
// dragging (>2 is too buggy for general use)
@@ -2441,7 +2445,7 @@ void GUI::create_option_window()
mesh_value[3]->align(FL_ALIGN_RIGHT);
mesh_value[3]->callback(mesh_options_ok_cb);
- mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 7 * BH, BW, BH, "Use incomplete high order elements (8-node quads, etc.)");
+ mesh_butt[4] = new Fl_Check_Button(L + 2 * WB, 2 * WB + 8 * BH, BW, BH, "Use incomplete high order elements (8-node quads, etc.)");
mesh_butt[4]->type(FL_TOGGLE_BUTTON);
mesh_butt[4]->callback(mesh_options_ok_cb);
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 0a45583fcceb29af4a9ae5c34d32ac866674253b..6cad537a008880179cba3c5c6ff48fbfd51bc585 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdge.cpp,v 1.30 2007-07-26 16:28:27 anand Exp $
+// $Id: GEdge.cpp,v 1.31 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -147,8 +147,6 @@ double GEdge::curvature(double par) const
GPoint P1 = point(par - eps1);
GPoint P2 = point(par + eps2);
-
-
double D = sqrt ((P1.x() - P2.x()) * (P1.x() - P2.x()) +
(P1.y() - P2.y()) * (P1.y() - P2.y()) +
(P1.z() - P2.z()) * (P1.z() - P2.z()));
diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 53d7388b95d639e6a6bdead9a6551e387829c8cd..eb51d58c221ced1759670d51ddace49a13971f44 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -1,4 +1,4 @@
-// $Id: GEdgeLoop.cpp,v 1.6 2007-09-04 13:47:01 remacle Exp $
+// $Id: GEdgeLoop.cpp,v 1.7 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -117,7 +117,7 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
GEdgeSigned *prevOne = 0;
- Msg(INFO,"Building a wire");
+ // Msg(INFO,"Building a wire");
GEdgeSigned ges(0,0);
while(wire.size()){
ges = nextOne(prevOne, wire);
@@ -126,7 +126,7 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
break;
}
prevOne = ⩾
- ges.print();
+ // ges.print();
loop.push_back(ges);
}
}
diff --git a/Geo/GEntity.h b/Geo/GEntity.h
index b4df6a9c599c54ec14a40749a15c22e3c0acad08..68d66bd808ef2abbc7a004b3e30606c2d6465734 100644
--- a/Geo/GEntity.h
+++ b/Geo/GEntity.h
@@ -72,6 +72,11 @@ class GEntity {
Line,
Circle,
Ellipse,
+ Conic,
+ Parabola,
+ Hyperbola,
+ TrimmedCurve,
+ OffsetCurve,
BSpline,
Bezier,
ParametricCurve,
@@ -105,6 +110,11 @@ class GEntity {
"Line",
"Circle",
"Ellipse",
+ "Conic",
+ "Parabola",
+ "Hyperbola",
+ "TrimmedCurve",
+ "OffsetCurve",
"BSpline",
"Bezier",
"Parametric curve",
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 3e422fd5b56cb8e9057b426f1e24b7565b4a3901..3865e9d9354c541f7b8d647ef35fb3c1874ef037 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.50 2007-12-15 03:41:57 geuzaine Exp $
+// $Id: GModel.cpp,v 1.51 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -529,7 +529,8 @@ static int checkVertices(std::vector<MVertex*> &vertices,
}
template <class T>
-static int checkElements(std::vector<T*> &elements,
+static int checkElements(int tag,
+ std::vector<T*> &elements,
std::set<MElement*, MElementLessThanLexicographic> &pos)
{
int num = 0;
@@ -540,8 +541,22 @@ static int checkElements(std::vector<T*> &elements,
pos.insert(e);
}
else{
- Msg(INFO, "Elements %d and %d have identical barycenter",
- (*it)->getNum(), e->getNum());
+ char temp[256];
+ char temp2[256];
+ sprintf(temp,"Elements %d tag %d(",(*it)->getNum(),tag);
+ for (int i=0;i<(*it)->getNumVertices();i++){
+ sprintf(temp2,"%d ",(*it)->getVertex(i)->getNum());
+ strcat(temp,temp2);
+ }
+ sprintf(temp2,") and %d(",e->getNum());
+ strcat(temp,temp2);
+ for (int i=0;i<e->getNumVertices();i++){
+ sprintf(temp2,"%d ",e->getVertex(i)->getNum());
+ strcat(temp,temp2);
+ }
+ sprintf(temp2,")have identical barycenter");
+ strcat(temp,temp2);
+ Msg(INFO, "%s",temp);
num++;
}
}
@@ -584,16 +599,16 @@ void GModel::checkMeshCoherence()
std::set<MElement*, MElementLessThanLexicographic> pos;
int num = 0;
for(eiter it = firstEdge(); it != lastEdge(); ++it)
- num += checkElements((*it)->lines, pos);
+ num += checkElements((*it)->tag(),(*it)->lines, pos);
for(fiter it = firstFace(); it != lastFace(); ++it){
- num += checkElements((*it)->triangles, pos);
- num += checkElements((*it)->quadrangles, pos);
+ num += checkElements((*it)->tag(),(*it)->triangles, pos);
+ num += checkElements((*it)->tag(),(*it)->quadrangles, pos);
}
for(riter it = firstRegion(); it != lastRegion(); ++it){
- num += checkElements((*it)->tetrahedra, pos);
- num += checkElements((*it)->hexahedra, pos);
- num += checkElements((*it)->prisms, pos);
- num += checkElements((*it)->pyramids, pos);
+ num += checkElements((*it)->tag(),(*it)->tetrahedra, pos);
+ num += checkElements((*it)->tag(),(*it)->hexahedra, pos);
+ num += checkElements((*it)->tag(),(*it)->prisms, pos);
+ num += checkElements((*it)->tag(),(*it)->pyramids, pos);
}
if(num) Msg(WARNING, "%d duplicate elements", num);
MElementLessThanLexicographic::tolerance = old_tol;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index a698a8a34cae90d28b06c03448d5b01d4c71f648..d0bcc47d5ad9059c56e39182b58a78864cdd4907 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -188,6 +188,7 @@ class GModel
int readOCCBREP(const std::string &name);
int readOCCIGES(const std::string &name);
int readOCCSTEP(const std::string &name);
+ void snapGVertices (void);
// Mesh IO
// =========================================
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 02638bef6203e9420256e6896d06d7f57b1f2069..b24ffdbac7f419e90cfd75fb79b3e245d9080710 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_OCC.cpp,v 1.22 2007-10-14 09:51:17 geuzaine Exp $
+// $Id: GModelIO_OCC.cpp,v 1.23 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -455,9 +455,42 @@ int GModel::readOCCBREP(const std::string &fn)
occ_internals->loadBREP(fn.c_str());
occ_internals->buildLists();
occ_internals->buildGModel(this);
+ snapGVertices();
return 1;
}
+void GModel::snapGVertices (void)
+{
+ viter vit = firstVertex();
+
+ double tol = CTX.geom.tolerance;
+
+ while (vit != lastVertex()){
+ std::list<GEdge*> edges = (*vit)->edges();
+ for (std::list<GEdge*>::iterator it = edges.begin() ;
+ it != edges .end(); ++it){
+ Range<double> parb = (*it)->parBounds(0);
+ double t;
+ if ((*it)->getBeginVertex() == *vit){
+ t = parb.low();
+ }
+ else if ((*it)->getEndVertex() == *vit){
+ t = parb.high();
+ }
+ else throw;
+ GPoint gp = (*it)->point(t);
+ double d = sqrt ((gp.x()-(*vit)->x())*(gp.x()-(*vit)->x())+
+ (gp.y()-(*vit)->y())*(gp.y()-(*vit)->y())+
+ (gp.z()-(*vit)->z())*(gp.z()-(*vit)->z()));
+ if (d > tol){
+ (*vit)->setPosition(gp);
+ Msg(WARNING, "Geom Vertex %d Corrupted (%12.5E)... Snap performed",(*vit)->tag(),d);
+ }
+ }
+ vit++;
+ }
+}
+
void GModel::deleteOCCInternals()
{
if(occ_internals) delete occ_internals;
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 25d1f22f99b0012476eef32a171491f0931b1df2..c2f072d2cef9cc1cb967f1e9b0173cf575ab071d 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoInterpolation.cpp,v 1.29 2007-09-03 20:09:14 geuzaine Exp $
+// $Id: GeoInterpolation.cpp,v 1.30 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -238,13 +238,14 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
V.u = u;
if(derivee) {
- double eps = 1.e-5;
+ double eps1 = u==0?0:1.e-5;
+ double eps2 = u==1?0:1.e-5;
Vertex D[2];
- D[0] = InterpolateCurve(c, u, 0);
- D[1] = InterpolateCurve(c, u + eps, 0);
- V.Pos.X = (D[1].Pos.X - D[0].Pos.X) / eps;
- V.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / eps;
- V.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / eps;
+ D[0] = InterpolateCurve(c, u-eps1, 0);
+ D[1] = InterpolateCurve(c, u+eps2, 0);
+ V.Pos.X = (D[1].Pos.X - D[0].Pos.X) / (eps1+eps2);
+ V.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / (eps1+eps2);
+ V.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / (eps1+eps2);
return V;
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 682897fc2a0b83eae82960416ba58ea4a20ae912..b8dd27abaffea06a324f035b37bccaa829c364f3 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.46 2007-11-26 14:34:09 remacle Exp $
+// $Id: MElement.cpp,v 1.47 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -687,14 +687,35 @@ double coef5[15][15]={
60.41666667, -38.54166667, 6.25000000}
};
+void GeomShapeFunctionP1(double u, double v, double *sf)
+{
+ for (int i = 0; i < 3; i++){
+ sf[i] = 0;
+ for(int j = 0; j < 3; j++){
+ sf[i] += coef1[i][j] * pow(u,P1[j][0]) * pow(v, P1[j][1]);
+ }
+ }
+}
+
+
void GradGeomShapeFunctionP1(double u, double v, double grads[6][2])
{
for (int i = 0; i < 3; i++){
grads[i][0] = 0;
grads[i][1] = 0;
for(int j = 0; j < 3; j++){
- if(P1[j][0] > 0) grads[i][0] += coef1[i][j] * pow(u,P1[j][0] - 1) * pow(v, P1[j][1]);
- if(P1[j][1] > 0) grads[i][1] += coef1[i][j] * pow(u,P1[j][0]) * pow(v, P1[j][1] - 1);
+ if(P1[j][0] > 0) grads[i][0] += coef1[i][j] * P1[j][0] * pow(u,P1[j][0] - 1) * pow(v, P1[j][1]);
+ if(P1[j][1] > 0) grads[i][1] += coef1[i][j] * P1[j][1] * pow(u,P1[j][0]) * pow(v, P1[j][1] - 1);
+ }
+ }
+}
+
+void GeomShapeFunctionP2(double u, double v, double *sf)
+{
+ for (int i = 0; i < 6; i++){
+ sf[i] = 0;
+ for(int j = 0; j < 6; j++){
+ sf[i] += coef2[i][j] * pow(u,P2[j][0]) * pow(v, P2[j][1]);
}
}
}
@@ -705,8 +726,18 @@ void GradGeomShapeFunctionP2(double u, double v, double grads[6][2])
grads[i][0] = 0;
grads[i][1] = 0;
for (int j = 0; j < 6; j++){
- if(P2[j][0] > 0) grads[i][0] += coef2[i][j] * pow(u, P2[j][0] - 1) * pow(v, P2[j][1]);
- if(P2[j][1] > 0) grads[i][1] += coef2[i][j] * pow(u, P2[j][0]) * pow(v, P2[j][1] - 1);
+ if(P2[j][0] > 0) grads[i][0] += coef2[i][j] * P2[j][0] * pow(u, P2[j][0] - 1) * pow(v, P2[j][1]);
+ if(P2[j][1] > 0) grads[i][1] += coef2[i][j] * P2[j][1] * pow(u, P2[j][0]) * pow(v, P2[j][1] - 1);
+ }
+ }
+}
+
+void GeomShapeFunctionP3(double u, double v, double *sf)
+{
+ for (int i = 0; i < 9; i++){
+ sf[i] = 0;
+ for(int j = 0; j < 9; j++){
+ sf[i] += coef3[i][j] * pow(u,P3[j][0]) * pow(v, P3[j][1]);
}
}
}
@@ -717,20 +748,41 @@ void GradGeomShapeFunctionP3 (double u, double v, double grads[9][2])
grads[i][0] = 0;
grads[i][1] = 0;
for(int j = 0; j < 9; j++){
- if(P3[j][0] > 0) grads[i][0] += coef3[i][j] * pow(u, P3[j][0] - 1) * pow(v, P3[j][1]);
- if(P3[j][1] > 0) grads[i][1] += coef3[i][j] * pow(u, P3[j][0]) * pow(v, P3[j][1] - 1);
+ if(P3[j][0] > 0) grads[i][0] += coef3[i][j] * P3[j][0] * pow(u, P3[j][0] - 1) * pow(v, P3[j][1]);
+ if(P3[j][1] > 0) grads[i][1] += coef3[i][j] * P3[j][1] * pow(u, P3[j][0]) * pow(v, P3[j][1] - 1);
}
}
}
+void GeomShapeFunctionP4(double u, double v, double *sf)
+{
+ for (int i = 0; i < 12; i++){
+ sf[i] = 0;
+ for(int j = 0; j < 12; j++){
+ sf[i] += coef4[i][j] * pow(u,P4[j][0]) * pow(v, P4[j][1]);
+ }
+ }
+}
+
+
void GradGeomShapeFunctionP4(double u, double v, double grads[12][2])
{
for(int i = 0; i < 12; i++){
grads[i][0] = 0;
grads[i][1] = 0;
for(int j = 0; j < 12; j++){
- if(P4[j][0] > 0) grads[i][0] += coef4[i][j] * pow(u, P4[j][0] - 1) * pow(v, P4[j][1]);
- if(P4[j][1] > 0) grads[i][1] += coef4[i][j] * pow(u, P4[j][0]) * pow(v, P4[j][1] - 1);
+ if(P4[j][0] > 0) grads[i][0] += coef4[i][j] * P4[j][0] * pow(u, P4[j][0] - 1) * pow(v, P4[j][1]);
+ if(P4[j][1] > 0) grads[i][1] += coef4[i][j] * P4[j][1] * pow(u, P4[j][0]) * pow(v, P4[j][1] - 1);
+ }
+ }
+}
+
+void GeomShapeFunctionP5(double u, double v, double *sf)
+{
+ for (int i = 0; i < 15; i++){
+ sf[i] = 0;
+ for(int j = 0; j < 15; j++){
+ sf[i] += coef5[i][j] * pow(u,P5[j][0]) * pow(v, P5[j][1]);
}
}
}
@@ -741,15 +793,17 @@ void GradGeomShapeFunctionP5(double u, double v, double grads[15][2])
grads[i][0] = 0;
grads[i][1] = 0;
for (int j = 0; j < 15; j++){
- if(P5[j][0] > 0) grads[i][0] += coef5[i][j] * pow(u, P5[j][0] - 1) * pow(v, P5[j][1]);
- if(P5[j][1] > 0) grads[i][1] += coef5[i][j] * pow(u, P5[j][0]) * pow(v, P5[j][1] - 1);
+ if(P5[j][0] > 0) grads[i][0] += coef5[i][j] * P5[j][0] * pow(u, P5[j][0] - 1) * pow(v, P5[j][1]);
+ if(P5[j][1] > 0) grads[i][1] += coef5[i][j] * P5[j][1] * pow(u, P5[j][0]) * pow(v, P5[j][1] - 1);
}
}
}
-void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double j[2][2])
+void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double j[2][3])
{
double grads[256][2];
+
+
switch(ord){
case 1: GradGeomShapeFunctionP1(uu, vv, grads); break;
case 2: GradGeomShapeFunctionP2(uu, vv, grads); break;
@@ -762,23 +816,135 @@ void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double j[2][2]
j[1][0] = 0 ; for(int i = 0; i < 3; i++) j[1][0] += grads [i][1] * _v[i] -> x();
j[0][1] = 0 ; for(int i = 0; i < 3; i++) j[0][1] += grads [i][0] * _v[i] -> y();
j[1][1] = 0 ; for(int i = 0; i < 3; i++) j[1][1] += grads [i][1] * _v[i] -> y();
+ j[0][2] = 0 ; for(int i = 0; i < 3; i++) j[0][2] += grads [i][0] * _v[i] -> z();
+ j[1][2] = 0 ; for(int i = 0; i < 3; i++) j[1][2] += grads [i][1] * _v[i] -> z();
+
+
for(int i = 3; i < 3 * ord; i++) j[0][0] += grads[i][0] * vs[i - 3] -> x();
for(int i = 3; i < 3 * ord; i++) j[1][0] += grads[i][1] * vs[i - 3] -> x();
for(int i = 3; i < 3 * ord; i++) j[0][1] += grads[i][0] * vs[i - 3] -> y();
for(int i = 3; i < 3 * ord; i++) j[1][1] += grads[i][1] * vs[i - 3] -> y();
+ for(int i = 3; i < 3 * ord; i++) j[0][2] += grads[i][0] * vs[i - 3] -> z();
+ for(int i = 3; i < 3 * ord; i++) j[1][2] += grads[i][1] * vs[i - 3] -> z();
+
+}
+
+void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, SPoint3 &p)
+{
+ double sf[256];
+
+ switch(ord){
+ case 1: GeomShapeFunctionP1(uu, vv, sf); break;
+ case 2: GeomShapeFunctionP2(uu, vv, sf); break;
+ case 3: GeomShapeFunctionP3(uu, vv, sf); break;
+ case 4: GeomShapeFunctionP4(uu, vv, sf); break;
+ case 5: GeomShapeFunctionP5(uu, vv, sf); break;
+ default: throw;
+ }
+
+ double x = 0 ; for(int i = 0; i < 3; i++) x += sf[i] * _v[i] -> x();
+ double y = 0 ; for(int i = 0; i < 3; i++) y += sf[i] * _v[i] -> y();
+ double z = 0 ; for(int i = 0; i < 3; i++) z += sf[i] * _v[i] -> z();
+
+ for(int i = 3; i < 3 * ord; i++) x += sf[i] * vs[i - 3] -> x();
+ for(int i = 3; i < 3 * ord; i++) y += sf[i] * vs[i - 3] -> y();
+ for(int i = 3; i < 3 * ord; i++) z += sf[i] * vs[i - 3] -> z();
+
+ p = SPoint3(x,y,z);
+
}
-void MTriangleN::jac(double uu, double vv , double j[2][2])
+
+
+void MTriangleN::jac(double uu, double vv , double j[2][3])
{
MTriangle::jac(_order, &(*(_vs.begin())), uu, vv, j);
}
-void MTriangle6::jac(double uu, double vv , double j[2][2])
+void MTriangleN::pnt(double uu, double vv, SPoint3 &p){
+ MTriangle::pnt(_order, &(*(_vs.begin())), uu, vv, p);
+}
+
+void MTriangle6::jac(double uu, double vv , double j[2][3])
{
MTriangle::jac(2, _vs, uu, vv, j);
}
-void MTriangle::jac(double uu, double vv, double j[2][2])
+void MTriangle6::pnt(double uu, double vv, SPoint3 &p){
+ MTriangle::pnt(2, _vs, uu, vv, p);
+}
+
+void MTriangle::jac(double uu, double vv, double j[2][3])
{
jac(1, 0, uu, vv, j);
}
+
+void MTriangle::pnt(double uu, double vv, SPoint3 &p){
+ MTriangle::pnt(1, 0, uu, vv, p);
+}
+
+int MTriangle6::getNumEdgesRep(){ return 30; }
+void MTriangle6::getEdgeRep (int num, double *x, double *y, double *z, SVector3 *n){
+ if (num < 10){
+ SPoint3 pnt1,pnt2;
+ pnt ( (double)num/10. , 0. , pnt1);
+ pnt ( (double)(num+1)/10. , 0. , pnt2);
+ x[0] = pnt1.x();x[1] = pnt2.x();
+ y[0] = pnt1.y();y[1] = pnt2.y();
+ z[0] = pnt1.z();z[1] = pnt2.z();
+ return;
+ }
+
+ if (num < 20){
+ SPoint3 pnt1,pnt2;
+ num -=10;
+ pnt ( 1.-(double)num/10. , (double)num/10. , pnt1);
+ pnt ( 1.-(double)(num+1)/10. , (double)(num+1)/10. , pnt2);
+ x[0] = pnt1.x();x[1] = pnt2.x();
+ y[0] = pnt1.y();y[1] = pnt2.y();
+ z[0] = pnt1.z();z[1] = pnt2.z();
+ return ;
+ }
+ {
+ SPoint3 pnt1,pnt2;
+ num -= 20;
+ pnt ( 0,(double)num/10. , pnt1);
+ pnt ( 0,(double)(num+1)/10., pnt2);
+ x[0] = pnt1.x();x[1] = pnt2.x();
+ y[0] = pnt1.y();y[1] = pnt2.y();
+ z[0] = pnt1.z();z[1] = pnt2.z();
+ }
+}
+
+int MTriangleN::getNumEdgesRep(){ return 120; }
+void MTriangleN::getEdgeRep (int num, double *x, double *y, double *z, SVector3 *n){
+ if (num < 40){
+ SPoint3 pnt1,pnt2;
+ pnt ( (double)num/40. , 0. , pnt1);
+ pnt ( (double)(num+1)/40. , 0. , pnt2);
+ x[0] = pnt1.x();x[1] = pnt2.x();
+ y[0] = pnt1.y();y[1] = pnt2.y();
+ z[0] = pnt1.z();z[1] = pnt2.z();
+ return;
+ }
+
+ if (num < 80){
+ SPoint3 pnt1,pnt2;
+ num -=40;
+ pnt ( 1.-(double)num/40. , (double)num/40. , pnt1);
+ pnt ( 1.-(double)(num+1)/40. , (double)(num+1)/40. , pnt2);
+ x[0] = pnt1.x();x[1] = pnt2.x();
+ y[0] = pnt1.y();y[1] = pnt2.y();
+ z[0] = pnt1.z();z[1] = pnt2.z();
+ return ;
+ }
+ {
+ SPoint3 pnt1,pnt2;
+ num -= 80;
+ pnt ( 0,(double)num/40. , pnt1);
+ pnt ( 0,(double)(num+1)/40., pnt2);
+ x[0] = pnt1.x();x[1] = pnt2.x();
+ y[0] = pnt1.y();y[1] = pnt2.y();
+ z[0] = pnt1.z();z[1] = pnt2.z();
+ }
+}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index f5d3fcbce7e19023fb2eb2ca5de6946adbfeab5a..f9c93b1436f0006c6be328215c96f3ec59c49bd0 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -304,8 +304,9 @@ class MLineN : public MLine {
class MTriangle : public MElement {
protected:
MVertex *_v[3];
- virtual void jac(int order, MVertex *verts[], double u, double v, double jac[2][2]);
public :
+ virtual void jac(int order, MVertex *verts[], double u, double v, double jac[2][3]);
+ virtual void pnt(int order, MVertex *verts[], double u, double v, SPoint3 &);
MTriangle(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0)
: MElement(num, part)
{
@@ -378,7 +379,8 @@ class MTriangle : public MElement {
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
}
- virtual void jac(double u, double v, double j[2][2]);
+ virtual void jac(double u, double v, double j[2][3]);
+ virtual void pnt(double u, double v, SPoint3&);
};
class MTriangle6 : public MTriangle {
@@ -408,16 +410,16 @@ class MTriangle6 : public MTriangle {
return getVertex(map[num]);
}
virtual int getNumEdgeVertices(){ return 3; }
- virtual int getNumEdgesRep(){ return 6; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[6][2] = {
- {0, 3}, {3, 1},
- {1, 4}, {4, 2},
- {2, 5}, {5, 0}
- };
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0);
- }
+ virtual int getNumEdgesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
+/* { */
+/* static const int e[6][2] = { */
+/* {0, 3}, {3, 1}, */
+/* {1, 4}, {4, 2}, */
+/* {2, 5}, {5, 0} */
+/* }; */
+/* _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, 0); */
+/* } */
virtual int getNumFacesRep(){ return 4; }
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -437,7 +439,8 @@ class MTriangle6 : public MTriangle {
tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
tmp = _vs[0]; _vs[0] = _vs[2]; _vs[2] = tmp;
}
- virtual void jac(double u, double v, double j[2][2]);
+ virtual void jac ( double u, double v , double j[2][3]) ;
+ virtual void pnt(double u, double v, SPoint3&);
};
class MTriangleN : public MTriangle {
@@ -481,13 +484,8 @@ class MTriangleN : public MTriangle {
throw;
}
virtual int getNumEdgeVertices(){ return _order - 1; }
- virtual int getNumEdgesRep(){ return 3 * _order; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- _getEdgeRep(getVertex(_orderedIndex(num)),
- getVertex(_orderedIndex((num + 1) % (3 * _order))),
- x, y, z, n, 0);
- }
+ virtual int getNumEdgesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual int getNumFacesRep(){ return 0; }
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -514,7 +512,8 @@ class MTriangleN : public MTriangle {
inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
_vs = inv;
}
- virtual void jac(double u, double v, double jac[2][2]);
+ virtual void jac(double u, double v, double jac[2][3]);
+ virtual void pnt(double u, double v, SPoint3&);
};
class MQuadrangle : public MElement {
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 9406a296fbc6affb6be383189a115adba1d6264e..c5c26189e39a4b1f07d7edee24530af3b266787f 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -124,8 +124,8 @@ class MEdgeVertex : public MVertex{
{
}
virtual ~MEdgeVertex(){}
- virtual bool getParameter(int i, double &par) const{ par = _u; return true; }
- virtual bool setParameter(int i, double par){ _u = par; return true; }
+ virtual bool getParameter(int i, double &par) const{ if(i)throw; par = _u; return true; }
+ virtual bool setParameter(int i, double par){ if(i)throw; _u = par; return true; }
double getLc () const {return _lc;}
};
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 95223af5c84936f8e29b7a2f392da5d131c68c59..c60cac8465ab0289bd90e5a5abb0cb0044782567 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCEdge.cpp,v 1.28 2007-12-03 15:17:40 remacle Exp $
+// $Id: OCCEdge.cpp,v 1.29 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -29,10 +29,14 @@ extern Context_T CTX;
#if defined(HAVE_OCC)
#include "Geom2dLProp_CLProps2d.hxx"
#include "Geom_BezierCurve.hxx"
-#include "Geom_BezierCurve.hxx"
+#include "Geom_OffsetCurve.hxx"
#include "Geom_Ellipse.hxx"
+#include "Geom_Parabola.hxx"
+#include "Geom_Hyperbola.hxx"
+#include "Geom_TrimmedCurve.hxx"
#include "Geom_Circle.hxx"
#include "Geom_Line.hxx"
+#include "Geom_Conic.hxx"
OCCEdge::OCCEdge(GModel *model, TopoDS_Edge edge, int num, GVertex *v1, GVertex *v2)
: GEdge(model, num, v1, v2), c(edge), trimmed(0)
@@ -89,12 +93,12 @@ SPoint2 OCCEdge::reparamOnFace(GFace *face, double epar,int dir) const
const double dz = p1.z()-p2.z();
if(sqrt(dx*dx+dy*dy+dz*dz) > 1.e-4 * CTX.lc){
// return reparamOnFace(face, epar,-1);
-// Msg(WARNING, "Reparam on face partially failed for curve %d surface %d at point %g",
-// tag(), face->tag(), epar);
-// Msg(WARNING, "On the face %d local (%g %g) global (%g %g %g)",
-// face->tag(), u, v, p2.x(), p2.y(), p2.z());
-// Msg(WARNING, "On the edge %d local (%g) global (%g %g %g)",
-// tag(), epar, p1.x(), p1.y(), p1.z());
+ Msg(WARNING, "Reparam on face partially failed for curve %d surface %d at point %g",
+ tag(), face->tag(), epar);
+ Msg(WARNING, "On the face %d local (%g %g) global (%g %g %g)",
+ face->tag(), u, v, p2.x(), p2.y(), p2.z());
+ Msg(WARNING, "On the edge %d local (%g) global (%g %g %g)",
+ tag(), epar, p1.x(), p1.y(), p1.z());
// GPoint ppp = face->closestPoint(SPoint3(p1.x(), p1.y(), p1.z()));
// return SPoint2(ppp.u(), ppp.v());
}
@@ -118,15 +122,15 @@ int OCCEdge::isSeam(GFace *face) const
GPoint OCCEdge::point(double par) const
{
- if(!curve.IsNull()){
- gp_Pnt pnt = curve->Value (par);
- return GPoint(pnt.X(), pnt.Y(), pnt.Z());
- }
- else if(trimmed){
+ if(trimmed){
double u, v;
curve2d->Value(par).Coord(u, v);
return trimmed->point(u, v);
}
+ else if(!curve.IsNull()){
+ gp_Pnt pnt = curve->Value (par);
+ return GPoint(pnt.X(), pnt.Y(), pnt.Z());
+ }
else{
Msg(WARNING, "OCC Curve %d is neither a 3D curve not a trimmed curve", tag());
return GPoint(0, 0, 0);
@@ -167,10 +171,20 @@ GEntity::GeomType OCCEdge::geomType() const
return Line;
else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_Ellipse))
return Ellipse;
+ else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_Parabola))
+ return Parabola;
+ else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_Hyperbola))
+ return Hyperbola;
+ else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve))
+ return TrimmedCurve;
+ else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve))
+ return OffsetCurve;
else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_BSplineCurve))
return BSpline;
else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_BezierCurve))
return Bezier;
+ else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_Conic))
+ return Conic;
return Unknown;
}
else{
@@ -178,12 +192,22 @@ GEntity::GeomType OCCEdge::geomType() const
return Circle;
else if (curve->DynamicType() == STANDARD_TYPE(Geom_Line))
return Line;
+ else if (curve->DynamicType() == STANDARD_TYPE(Geom_Parabola))
+ return Parabola;
+ else if (curve->DynamicType() == STANDARD_TYPE(Geom_Hyperbola))
+ return Hyperbola;
+ else if (curve->DynamicType() == STANDARD_TYPE(Geom_TrimmedCurve))
+ return TrimmedCurve;
+ else if (curve->DynamicType() == STANDARD_TYPE(Geom_OffsetCurve))
+ return OffsetCurve;
else if (curve->DynamicType() == STANDARD_TYPE(Geom_Ellipse))
return Ellipse;
else if (curve->DynamicType() == STANDARD_TYPE(Geom_BSplineCurve))
return BSpline;
else if (curve->DynamicType() == STANDARD_TYPE(Geom_BezierCurve))
return Bezier;
+ else if (curve2d->DynamicType() == STANDARD_TYPE(Geom_Conic))
+ return Conic;
return Unknown;
}
}
diff --git a/Geo/OCCVertex.h b/Geo/OCCVertex.h
index cc4150467153c171ecb3c79ce46bc1274abfb124..49bd651935169a8dbbe8422f90181206918055e9 100644
--- a/Geo/OCCVertex.h
+++ b/Geo/OCCVertex.h
@@ -30,12 +30,17 @@
class OCCVertex : public GVertex {
protected:
TopoDS_Vertex v;
+ double _x,_y,_z;
mutable double max_curvature;
double max_curvature_of_surfaces() const;
public:
OCCVertex(GModel *m, int num, TopoDS_Vertex _v) : GVertex(m, num), v(_v)
{
max_curvature = -1;
+ gp_Pnt pnt = BRep_Tool::Pnt (v);
+ _x = pnt.X();
+ _y = pnt.Y();
+ _z = pnt.Z();
mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
}
virtual ~OCCVertex() {}
@@ -45,19 +50,25 @@ class OCCVertex : public GVertex {
}
virtual double x() const
{
- gp_Pnt pnt = BRep_Tool::Pnt (v);
- return pnt.X();
+ return _x;
}
virtual double y() const
{
- gp_Pnt pnt = BRep_Tool::Pnt (v);
- return pnt.Y();
+ return _y;
}
virtual double z() const
{
- gp_Pnt pnt = BRep_Tool::Pnt (v);
- return pnt.Z();
+ return _z;
}
+ virtual void setPosition(GPoint &p){
+ _x = p.x();
+ _y = p.y();
+ _z = p.z();
+ delete mesh_vertices[0];
+ mesh_vertices.clear();
+ mesh_vertices.push_back(new MVertex(x(), y(), z(), this));
+ }
+
ModelType getNativeType() const { return OpenCascadeModel; }
void * getNativePtr() const { return (void*) &v; }
virtual SPoint2 reparamOnFace ( GFace *gf , int) const;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index bb5ad7cf6c1f6862cecfaeeb2c94bd140e30b3c9..107586b149c268302cae6f165bedfcc999ceed82 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -1,4 +1,4 @@
-// $Id: BDS.cpp,v 1.86 2007-11-26 14:34:09 remacle Exp $
+// $Id: BDS.cpp,v 1.87 2008-01-14 21:29:13 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -44,7 +44,7 @@ void outputScalarField(std::list < BDS_Face * >t, const char *iii, int param)
fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
pts[0]->u, pts[0]->v, 0.0,
pts[1]->u, pts[1]->v, 0.0,
- pts[2]->u, pts[2]->v, 0.0,(double)pts[0]->iD,(double)pts[1]->iD,(double)pts[2]->iD);
+ pts[2]->u, pts[2]->v, 0.0,(double)pts[0]->lc(),(double)pts[1]->lc(),(double)pts[2]->lc());
else
fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
@@ -1166,6 +1166,103 @@ bool test_move_point_parametric_triangle (BDS_Point * p, double u, double v, BDS
return ori_init*ori_final > 0;
}
+
+/**
+ d^2_i = (x^2_i - x)^T M (x_i - x)
+ = M11 (x_i - x)^2 + 2 M21 (x_i-x)(y_i-y) + M22 (y_i-y)^2
+
+ ~:-)
+
+
+*/
+
+struct smoothVertexData{
+ BDS_Point *p;
+ GFace *gf;
+ double scalu,scalv;
+ std::list < BDS_Face * >ts;
+};
+
+double smoothing_objective_function(double U, double V, BDS_Point *v, std::list < BDS_Face * >&ts,double su, double sv,GFace *gf){
+
+ GPoint gp = gf->point(U*su,V*sv);
+
+ const double oldX = v->X;
+ const double oldY = v->Y;
+ const double oldZ = v->Z;
+ v->X = gp.x();
+ v->Y = gp.y();
+ v->Z = gp.z();
+
+ std::list < BDS_Face * >::iterator it = ts.begin();
+ std::list < BDS_Face * >::iterator ite = ts.end();
+ double qMin=1;
+ while(it != ite) {
+ BDS_Face *t = *it;
+ qMin = std::min(qmTriangle(*it,QMTRI_RHO),qMin);
+ ++it;
+ }
+ v->X = oldX;
+ v->Y = oldY;
+ v->Z = oldZ;
+ return -qMin;
+}
+
+void deriv_smoothing_objective_function(double U, double V,
+ double &F, double &dFdU, double &dFdV,void *data){
+ smoothVertexData *svd = (smoothVertexData*)data;
+ BDS_Point *v = svd->p;
+ const double LARGE = 1.e5;
+ const double SMALL = 1./LARGE;
+ F = smoothing_objective_function(U,V,v,svd->ts,svd->scalu,svd->scalv,svd->gf);
+ double F_U = smoothing_objective_function(U+SMALL,V,v,svd->ts,svd->scalu,svd->scalv,svd->gf);
+ double F_V = smoothing_objective_function(U,V+SMALL,v,svd->ts,svd->scalu,svd->scalv,svd->gf);
+ dFdU = (F_U-F)*LARGE;
+ dFdV = (F_V-F)*LARGE;
+}
+
+double smooth_obj(double U, double V, void *data){
+ smoothVertexData *svd = (smoothVertexData*)data;
+ return smoothing_objective_function(U,V,svd->p,svd->ts,svd->scalu,svd->scalv,svd->gf);
+}
+
+
+void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv){
+#ifdef HAVE_GSL
+ if(data->g && data->g->classif_degree <= 1)
+ return;
+ smoothVertexData vd;
+ vd.p = data;
+ vd.scalu = su;
+ vd.scalv = sv;
+ vd.gf = GF;
+ data->getTriangles(vd.ts);
+ double U=data->u,V=data->v,val;
+
+ val = smooth_obj(U, V, &vd);
+ if (val < -.90)return;
+
+ minimize_2 ( smooth_obj, deriv_smoothing_objective_function, &vd, 5, U,V,val);
+ std::list < BDS_Face * >::iterator it = vd.ts.begin();
+ std::list < BDS_Face * >::iterator ite = vd.ts.end();
+ while(it != ite) {
+ BDS_Face *t = *it;
+ if (!test_move_point_parametric_triangle ( data, U, V, t)){
+ return;
+ }
+ ++it;
+ }
+
+ data->u = U;
+ data->v = V;
+ GPoint gp = GF->point(U*su,V*sv);
+ data->X = gp.x();
+ data->Y = gp.y();
+ data->Z = gp.z();
+#endif
+}
+
+
bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf)
{
@@ -1193,8 +1290,8 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf)
BDS_Face *t = *it;
BDS_Point *n[4];
t->getNodes(n);
- double S = fabs(surface_triangle(n[0],n[1],n[2]));
- S = 1;
+ // double S = fabs(surface_triangle(n[0],n[1],n[2]));
+ double S = 1;
sTot += S;
U += (n[0]->u + n[1]->u + n[2]->u) *S;
V += (n[0]->v + n[1]->v + n[2]->v) *S;
@@ -1211,14 +1308,25 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf)
const double oldY = p->Y;
const double oldZ = p->Z;
- double oldWorst=1.;
- double newWorst=1.;
+ double oldU=U;
+ double oldV=V;
it = ts.begin();
+ double s1=0,s2=0;
while(it != ite) {
BDS_Face *t = *it;
- if (!test_move_point_parametric_triangle ( p, U, V, t))
- return false;
+ BDS_Point *n[4];
+ t->getNodes(n);
+ p->u = U;
+ p->v = V;
+ s1 += fabs(surface_triangle_param(n[0],n[1],n[2]));
+ p->u = oldU;
+ p->v = oldV;
+ s2 += fabs(surface_triangle_param(n[0],n[1],n[2]));
+
+// if (!test_move_point_parametric_triangle ( p, U, V, t)){
+// return false;
+// }
// p->X = gp.x();
// p->Y = gp.y();
// p->Z = gp.z();
@@ -1230,6 +1338,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point * p, GFace *gf)
++it;
}
+ if (fabs(s2-s1) > 1.e-10 * (s2+s1))return false;
// if (newWorst < 1.e-2)
// {
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index bc7aa40110c8c0cebb4ac351f8280f562e7a93cf..9019311c0bb14ed30167964e824bddfdc8dea096 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -45,6 +45,7 @@ void vector_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3);
+void optimize_vertex_position (GFace *GF, BDS_Point *data, double su, double sv);
class BDS_GeomEntity
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 08a2dcb922e227b1f37cd9f9e543dc1bfb55725e..1a8f280bb3a8e47e006541e619e21db27dabf3c2 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -1,4 +1,4 @@
-// $Id: BackgroundMesh.cpp,v 1.30 2007-11-26 14:34:09 remacle Exp $
+// $Id: BackgroundMesh.cpp,v 1.31 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -70,25 +70,6 @@ static double max_edge_curvature(const GVertex *gv)
return max_curvature;
}
-static double max_surf_curvature(const GVertex *gv)
-{
- double max_curvature = 0;
- std::list<GEdge*> l_edges = gv->edges();
- for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
- ite != l_edges.end(); ++ite){
- GEdge *_myGEdge = *ite;
- std::list<GFace *> faces = _myGEdge->faces();
- std::list<GFace *>::iterator it = faces.begin();
- while (it != faces.end()){
- SPoint2 par = gv->reparamOnFace((*it),1);
- double cc = (*it)->curvature ( par );
- max_curvature = std::max ( cc , max_curvature);
- ++it;
- }
- }
- return max_curvature;
-}
-
static double max_surf_curvature(const GEdge *ge, double u)
{
double max_curvature = 0;
@@ -103,6 +84,23 @@ static double max_surf_curvature(const GEdge *ge, double u)
return max_curvature;
}
+static double max_surf_curvature(const GVertex *gv)
+{
+ double max_curvature = 0;
+ std::list<GEdge*> l_edges = gv->edges();
+ for (std::list<GEdge*>::const_iterator ite = l_edges.begin();
+ ite != l_edges.end(); ++ite){
+ GEdge *_myGEdge = *ite;
+ Range<double> bounds = _myGEdge->parBounds(0);
+ if (gv == _myGEdge->getBeginVertex())
+ max_curvature = std::max(max_curvature, max_surf_curvature(_myGEdge,bounds.low()));
+ else
+ max_curvature = std::max(max_curvature, max_surf_curvature(_myGEdge,bounds.high()));
+ }
+ return max_curvature;
+}
+
+
// the mesh vertex is classified on a model vertex. we compute the
// maximum of the curvature of model faces surrounding this point if
// it is classified on a model edge, we do the same for all model
@@ -114,17 +112,18 @@ double LC_MVertex_CURV(GEntity *ge, double U, double V)
double Crv = 0;
switch(ge->dim()){
case 0:
- Crv = max_edge_curvature ( (const GVertex *)ge);
- Crv = std::max(max_surf_curvature ( (const GVertex *)ge),Crv);
+ // Crv = max_edge_curvature ( (const GVertex *)ge);
+ // Crv = std::max(max_surf_curvature ( (const GVertex *)ge),Crv);
+ Crv = max_surf_curvature ( (const GVertex *)ge);
// printf("point %d coucou %g\n",ge->tag(),Crv);
break;
case 1:
{
GEdge *ged = (GEdge *)ge;
- Crv = ged->curvature(U);
+ // Crv = ged->curvature(U);
// printf("coucou %12.5E %d\n",Crv,CTX.mesh.min_circ_points);
- Crv = std::max(Crv,max_surf_curvature(ged, U));
-
+ // Crv = std::max(Crv,max_surf_curvature(ged, U));
+ Crv = max_surf_curvature(ged, U);
}
break;
case 2:
@@ -194,15 +193,16 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
lc = std::min(std::min(std::min(l1, l2), l3), l4) * CTX.mesh.lc_factor;
}
-
- if(CTX.mesh.lc_from_curvature && ge->dim() < 3)
+ if(CTX.mesh.lc_from_curvature && ge->dim() <= 2 )
lc = std::min (lc,LC_MVertex_CURV(ge, U, V));
+ lc = std::max(lc,CTX.mesh.lc_min*CTX.mesh.lc_factor);
+
if(lc <= 0.){
Msg(GERROR, "Incorrect char. length lc = %g: using default instead", lc);
return l3 * CTX.mesh.lc_factor;
}
-
+
return lc;
}
@@ -211,7 +211,7 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
// we do it also if CTX.mesh.constrained_bgmesh is true;
bool Extend1dMeshIn2dSurfaces()
{
- if(lc_field.empty() && !CTX.mesh.lc_from_curvature) return true;
+ if(lc_field.empty()) return true;
if(CTX.mesh.constrained_bgmesh) return true;
return false;
}
diff --git a/Mesh/BackgroundMesh.h b/Mesh/BackgroundMesh.h
index ec7958c8cb3a1842e086f935a6b8cf02e6969b25..30bd8facee1771732bc117472e47321b157a3161 100644
--- a/Mesh/BackgroundMesh.h
+++ b/Mesh/BackgroundMesh.h
@@ -28,4 +28,6 @@ void BGMAddField(Field *field);
void BGMReset();
bool Extend1dMeshIn2dSurfaces ();
bool Extend2dMeshIn3dVolumes ();
+double LC_MVertex_CURV(GEntity *ge, double U, double V);
+double LC_MVertex_PNTS(GEntity *ge, double U, double V);
#endif
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index 2c66a453c1c4f7ca33a89c4cd17ad575dea3a1ae..ab3d73428b3ffc07e9531df78f45a71cc60f08cc 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -1,4 +1,4 @@
-// $Id: Field.cpp,v 1.8 2007-11-09 08:07:53 geuzaine Exp $
+// $Id: Field.cpp,v 1.9 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -26,6 +26,7 @@
#include "Field.h"
#include "Context.h"
#include "GeoInterpolation.h"
+#include "BackgroundMesh.h"
#include "GModel.h"
#ifdef HAVE_MATH_EVAL
#include "matheval.h"
@@ -448,6 +449,50 @@ AttractorField_1DMesh::AttractorField_1DMesh(GModel *m, double dmax, double dmin
}
}
+AttractorField_1DMesh::AttractorField_1DMesh(GFace *gf, double dmax, double dmin,
+ double lcmax)
+ : _dmax(dmax), _dmin(dmin), _lcmax(lcmax)
+{
+ std::list<GEdge*> edges = gf->edges();
+ std::list<GEdge*>::iterator it = edges.begin();
+
+ std::map<MVertex*, double> maplc;
+ std::map<MVertex*, double> maplc2;
+
+ while(it != edges.end()){
+ MVertex *first = (*it)->getBeginVertex()->mesh_vertices[0];
+ for(unsigned int i = 1; i <= (*it)->mesh_vertices.size(); ++i){
+ MVertex *last = (i == (*it)->mesh_vertices.size()) ?
+ (*it)->getEndVertex()->mesh_vertices[0] : (*it)->mesh_vertices[i];
+ double l = sqrt((first->x() - last->x()) * (first->x() - last->x()) +
+ (first->y() - last->y()) * (first->y() - last->y()) +
+ (first->z() - last->z()) * (first->z() - last->z()));
+ double t;
+ last->getParameter(0,t);
+ double l2 = LC_MVertex_PNTS(last->onWhat(),t,0);
+ addMapLc(maplc, first, l);
+ addMapLc(maplc, last, l);
+ addMapLc(maplc2, first, l2);
+ addMapLc(maplc2, last, l2);
+ first = last;
+ }
+ ++it;
+ }
+
+ std::map<MVertex*, double>::iterator itm = maplc.begin();
+ while(itm != maplc.end()){
+ addPoint(itm->first->x(), itm->first->y(), itm->first->z());
+ lcs.push_back(itm->second);
+ ++itm;
+ }
+ itm = maplc2.begin();
+ while(itm != maplc2.end()){
+ lcs2.push_back(itm->second);
+ ++itm;
+ }
+}
+
+
double AttractorField_1DMesh::operator()(double X, double Y, double Z)
{
#ifdef HAVE_ANN_
@@ -463,3 +508,16 @@ double AttractorField_1DMesh::operator()(double X, double Y, double Z)
Msg(GERROR,"GMSH should be compiled with ANN in order to enable attractors");
#endif
}
+
+void AttractorField_1DMesh::eval(double X, double Y, double Z, double &lcmin, double &lcpt, double &d)
+{
+#ifdef HAVE_ANN_
+ double xyz[3] = {X, Y, Z};
+ kdtree->annkSearch(xyz, maxpts, index, dist);
+ d = sqrt(dist[0]);
+ lcmin = lcs[index[0]];
+ lcpt = lcs2[index[0]] * CTX.mesh.lc_factor;
+#else
+ Msg(GERROR,"GMSH should be compiled with ANN in order to enable attractors");
+#endif
+}
diff --git a/Mesh/Field.h b/Mesh/Field.h
index 9a0ce6bf35ad2f8327a2a47a4abe770d49d71991..723350ef94c11b00dd306014f12562901b734eba 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -139,10 +139,13 @@ class AttractorField_1DMesh : public AttractorField
{
protected:
std::vector<double> lcs;
+ std::vector<double> lcs2;
double _dmax,_dmin,_lcmax;
public:
- AttractorField_1DMesh (GModel *m, double dmax, double dmin, double lcmax);
+ AttractorField_1DMesh (GModel *m , double dmax, double dmin, double lcmax);
+ AttractorField_1DMesh (GFace *gf, double dmax, double dmin, double lcmax);
virtual double operator()(double X, double Y, double Z) ;
+ virtual void eval(double X, double Y, double Z, double &l, double &lpt, double &dist) ;
};
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index abb93ba35f4a6a41d5aae53721292d9b4168331e..13632b7017fa9836f384dae8726056a67f79f148 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.15 2007-09-12 20:14:34 geuzaine Exp $
+// $Id: HighOrder.cpp,v 1.16 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -20,6 +20,7 @@
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "HighOrder.h"
+#include "meshGFaceOptimize.h"
#include "MElement.h"
#include "Message.h"
#include "OS.h"
@@ -246,6 +247,7 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
uc < u0 || uc > u1){ // need to treat periodic curves properly!
SPoint3 pc = edge.interpolate(t);
v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
+ v->setParameter (0,t);
}
else {
GPoint pc = ge->point(uc);
@@ -379,6 +381,8 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
double vc = (1. - t1 - t2) * p0[1] + t1 * p1[1] + t2 * p2[1];
GPoint pc = gf->point(uc, vc);
v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
+ v->setParameter (0,uc);
+ v->setParameter (1,vc);
}
faceVertices[p].push_back(v);
gf->mesh_vertices.push_back(v);
@@ -698,28 +702,285 @@ bool straightLine(std::vector<MVertex*> &l, MVertex *n1, MVertex *n2)
double by = n1->y();
double ay = n2->y() - by;
double y = ay * t + by;
- if(fabs(y-v->y()) > 1.e-11 * CTX.lc){
+ if(fabs(y-v->y()) > 1.e-07 * CTX.lc){
return false;
}
}
return true;
}
+FILE *MYFILE = 0;
+
void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
{
- double j[2][2];
- int n = 3;
+ double mat[2][3];
+ int n = 10;
+
+ t->jac(1,0,0,0,mat);
+ double v1[3] = {mat[0][0],mat[0][1],mat[0][2]};
+ double v2[3] = {mat[1][0],mat[1][1],mat[1][2]};
+ double normal1[3],normal[3];
+ prodve(v1,v2,normal1);
+ norme(normal1);
+ double sign = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0];
+ double invsign = 1./(sign);
+
for(int i = 0; i < n; i++){
for(int k = 0; k < n - i; k++){
- t->jac((double)i / (n - 1), (double)k / (n - 1), j);
- double det = det2x2(j);
+ t->jac((double)i / (n - 1), (double)k / (n - 1), mat);
+ // printf("%g %g\n",(double)i / (n - 1), (double)k / (n - 1));
+ // SPoint3 pt;
+ // t->pnt((double)i / (n - 1), (double)k / (n - 1),pt);
+ double v1b[3] = {mat[0][0],mat[0][1],mat[0][2]};
+ double v2b[3] = {mat[1][0],mat[1][1],mat[1][2]};
+ prodve(v1b,v2b,normal);
+ //double sign; prosca(normal1,normal,&sign);
+ //double det = norm3(normal) * (sign>0?1.:-1.);
+ double det = invsign*(mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]) * (sign>0?1.:-1.);
+ // double det2 = 1./det1;
+ // double det = std::max(det1,det2);
+ // if(MYFILE)fprintf(MYFILE,"SP(%g,%g,%g){%g};\n",pt.x(),pt.y(),pt.z(),det);
minJ = std::min(det, minJ);
maxJ = std::max(det, maxJ);
}
}
}
-void smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
+struct smoothVertexDataHO{
+ MVertex *v;
+ GFace *gf;
+ std::vector < MTriangle * >ts;
+};
+
+struct smoothVertexDataHON{
+ std::vector<MVertex*> v;
+ GFace *gf;
+ std::vector < MTriangle * >ts;
+};
+
+double smoothing_objective_function_HighOrder (double U, double V, MVertex *v, std::vector < MTriangle * >&ts ,GFace *gf){
+
+ GPoint gp = gf->point(U,V);
+ const double oldX = v->x();
+ const double oldY = v->y();
+ const double oldZ = v->z();
+
+ v->x() = gp.x();
+ v->y() = gp.y();
+ v->z() = gp.z();
+
+ double minJ = 1.e22;
+ double maxJ = -1.e22;
+ for (int i=0;i<ts.size();i++){
+ getMinMaxJac (ts[i], minJ, maxJ);
+ }
+ v->x() = oldX;
+ v->y() = oldY;
+ v->z() = oldZ;
+
+ return -minJ;
+}
+
+void deriv_smoothing_objective_function_HighOrder(double U, double V,
+ double &F, double &dFdU, double &dFdV,void *data){
+ smoothVertexDataHO *svd = (smoothVertexDataHO*)data;
+ MVertex *v = svd->v;
+ const double LARGE = -1.e5;
+ const double SMALL = 1./LARGE;
+ F = smoothing_objective_function_HighOrder(U,V,v,svd->ts,svd->gf);
+ double F_U = smoothing_objective_function_HighOrder(U+SMALL,V,v,svd->ts,svd->gf);
+ double F_V = smoothing_objective_function_HighOrder(U,V+SMALL,v,svd->ts,svd->gf);
+ dFdU = (F_U-F)*LARGE;
+ dFdV = (F_V-F)*LARGE;
+}
+
+double smooth_obj_HighOrder(double U, double V, void *data){
+ smoothVertexDataHO *svd = (smoothVertexDataHO*)data;
+ return smoothing_objective_function_HighOrder(U,V,svd->v,svd->ts,svd->gf);
+}
+
+double smooth_obj_HighOrderN(double *uv, void *data){
+ smoothVertexDataHON *svd = (smoothVertexDataHON*)data;
+ double oldX[10],oldY[10],oldZ[10];
+ for (int i=0;i<svd->v.size();i++){
+ GPoint gp = svd->gf->point(uv[2*i],uv[2*i+1]);
+ oldX[i] = svd->v[i]->x();
+ oldY[i] = svd->v[i]->y();
+ oldZ[i] = svd->v[i]->z();
+ svd->v[i]->x() = gp.x();
+ svd->v[i]->y() = gp.y();
+ svd->v[i]->z() = gp.z();
+ }
+ double minJ = 1.e22;
+ double maxJ = -1.e22;
+ for (int i=0;i<svd->ts.size();i++){
+ getMinMaxJac (svd->ts[i], minJ, maxJ);
+ }
+ for (int i=0;i<svd->v.size();i++){
+ svd->v[i]->x() = oldX[i];
+ svd->v[i]->y() = oldY[i];
+ svd->v[i]->z() = oldZ[i];
+ }
+ return -minJ;
+}
+
+
+void deriv_smoothing_objective_function_HighOrderN(double *uv, double *dF, double &F,void *data){
+ const double LARGE = -1.e2;
+ const double SMALL = 1./LARGE;
+ smoothVertexDataHON *svd = (smoothVertexDataHON*)data;
+ F = smooth_obj_HighOrderN(uv,data);
+ for (int i=0;i<svd->v.size();i++){
+ uv[i] += SMALL;
+ dF[i] = (smooth_obj_HighOrderN(uv,data) - F)*LARGE;
+ uv[i] -= SMALL;
+ }
+}
+
+
+
+bool optimizeHighOrderMesh(GModel *modl, GFace *gf, edgeContainer &edgeVertices){
+
+ v2t_cont adjv;
+ buildVertexToTriangle ( gf->triangles , adjv );
+
+ int ITER=0;
+
+ typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MElement*> > edge2tris;
+ edge2tris e2t;
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ MTriangle *t = gf->triangles[i];
+ for(int j = 0; j < t->getNumEdges(); j++){
+ MEdge edge = t->getEdge(j);
+ std::pair<MVertex*, MVertex*> p(edge.getMinVertex(), edge.getMaxVertex());
+ e2t[p].push_back(t);
+ }
+ }
+
+ bool success = false;
+
+ for(edge2tris::iterator it = e2t.begin(); it != e2t.end(); ++it){
+ std::pair<MVertex*, MVertex*> edge = it->first;
+ std::vector<MVertex*> e;
+ std::vector<MElement*> triangles = it->second;
+ if(triangles.size() == 2){
+ MVertex *n2 = edge.first;
+ MVertex *n4 = edge.second;
+ MTriangle *t1 = (MTriangle*)triangles[0];
+ MTriangle *t2 = (MTriangle*)triangles[1];
+ if(n2 < n4)
+ e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n2, n4)];
+ else
+ e = edgeVertices[std::make_pair<MVertex*, MVertex*> (n4, n2)];
+
+ if (e.size() == -1){
+ double initu,initv;
+ e[0]->getParameter ( 0,initu);
+ e[0]->getParameter ( 1,initv);
+ smoothVertexDataHO vd;
+ vd.v = e[0];
+ vd.gf = gf;
+ vd.ts.push_back(t1);
+ vd.ts.push_back(t2);
+ double val;
+ minimize_2 ( smooth_obj_HighOrder, deriv_smoothing_objective_function_HighOrder, &vd, 1, initu,initv,val);
+ vd.v->setParameter(0,initu);
+ vd.v->setParameter(1,initv);
+ GPoint gp = gf->point(initu,initv);
+ vd.v->x() = gp.x();
+ vd.v->y() = gp.y();
+ vd.v->z() = gp.z();
+ }
+ else if (e.size() < 5){
+ double uv[20];
+ for (int i=0;i<e.size();i++){
+ if (!e[i]->getParameter (0,uv[2*i]))throw;
+ if (!e[i]->getParameter (1,uv[2*i+1]))throw;
+ }
+ smoothVertexDataHON vdN;
+ vdN.v = e;
+ vdN.gf = gf;
+ vdN.ts.clear();
+ vdN.ts.push_back(t1);
+ vdN.ts.push_back(t2);
+ double val;
+ double F = -smooth_obj_HighOrderN(uv, &vdN);
+ // printf("F = %12.5E %p %p\n",F,t1,t2);
+ if (F < .2){
+ minimize_N ( 2*e.size(), smooth_obj_HighOrderN, deriv_smoothing_objective_function_HighOrderN, &vdN, 1, uv,val);
+ double Fafter = -smooth_obj_HighOrderN(uv, &vdN);
+ if (F < Fafter){
+ // printf("found a pattern with f = %22.15E -> %22.15E (%22.15E) %d points\n",F,Fafter,val,e.size());
+ success = true;
+// checkHighOrderTriangles(modl);
+// double minJ=1.e22,maxJ=-1.e22;
+// getMinMaxJac (t1, minJ, maxJ);
+// printf("AVANT 1 minJ %22.15E maxJ %22.15E\n",minJ,maxJ);
+// minJ=1.e22;maxJ=-1.e22;
+// getMinMaxJac (t2, minJ, maxJ);
+// printf("AVANT 2 minJ %22.15E maxJ %22.15E\n",minJ,maxJ);
+// checkHighOrderTriangles(modl);
+
+ for (int i=0;i<vdN.v.size();i++){
+ vdN.v[i]->setParameter ( 0,uv[2*i]);
+ vdN.v[i]->setParameter ( 1,uv[2*i+1]);
+ GPoint gp = gf->point(uv[2*i],uv[2*i+1]);
+ vdN.v[i]->x() = gp.x();
+ vdN.v[i]->y() = gp.y();
+ vdN.v[i]->z() = gp.z();
+ }
+// minJ=1.e22;maxJ=-1.e22;
+// getMinMaxJac (t1, minJ, maxJ);
+// printf("APRES 1 minJ %22.15E maxJ %22.15E\n",minJ,maxJ);
+// minJ=1.e22;maxJ=-1.e22;
+// getMinMaxJac (t2, minJ, maxJ);
+// printf("APRES 2 minJ %22.15E maxJ %22.15E\n",minJ,maxJ);
+// checkHighOrderTriangles(modl);
+// return true;
+ }
+ }
+ }
+ }
+ }
+
+ while (1){
+ v2t_cont :: iterator it = adjv.begin();
+ while (it != adjv.end()){
+ MVertex *ver= it->first;
+ GEntity *ge = ver->onWhat();
+ if (ge->dim() == 2){
+ double initu,initv;
+ ver->getParameter ( 0,initu);
+ ver->getParameter ( 1,initv);
+ smoothVertexDataHO vd;
+ vd.v = ver;
+ vd.gf = gf;
+ vd.ts = it->second;
+ double val;
+
+ double F = -smooth_obj_HighOrder(initu,initv, &vd);
+ if (F < .2){
+ minimize_2 ( smooth_obj_HighOrder, deriv_smoothing_objective_function_HighOrder, &vd, 1, initu,initv,val);
+ double Fafter = -smooth_obj_HighOrder(initu,initv, &vd);
+ if (F < Fafter){
+ success = true;
+ ver->setParameter(0,initu);
+ ver->setParameter(1,initv);
+ GPoint gp = gf->point(initu,initv);
+ ver->x() = gp.x();
+ ver->y() = gp.y();
+ ver->z() = gp.z();
+ }
+ }
+ }
+ ++it;
+ }
+ break;
+ }
+ return success;
+}
+
+bool smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
{
typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MElement*> > edge2tris;
edge2tris e2t;
@@ -732,6 +993,10 @@ void smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
}
}
+ bool success = false;
+
+ const int NBST = 10;
+
for(edge2tris::iterator it = e2t.begin(); it != e2t.end(); ++it){
std::pair<MVertex*, MVertex*> edge = it->first;
std::vector<MVertex*> e1, e2, e3, e4, e;
@@ -766,68 +1031,123 @@ void smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
if((!straightLine(e1, n1, n2) || !straightLine(e2, n2, n3) ||
!straightLine(e3, n3, n4) || !straightLine(e4, n4, n1))){
+
+
+ double Unew[NBST][10],Vnew[NBST][10];
+ double Xold[10],Yold[10],Zold[10];
+
for(unsigned int i = 0; i < e.size(); i++){
- double v = (double)(i + 1) / (e.size() + 1);
- double u = 1. - v;
- MVertex *vert = (n2 < n4) ? e[i] : e[e.size() - i - 1];
- MVertex *vert1 = (n1 < n2) ? e1[e1.size() - i - 1] : e1[i];
- MVertex *vert3 = (n3 < n4) ? e3[i] : e3[e3.size() - i - 1];
- MVertex *vert4 = (n4 < n1) ? e4[e4.size() - i - 1] : e4[i];
- MVertex *vert2 = (n2 < n3) ? e2[i] : e2[e2.size() - i - 1];
- vert->x() = vert->x() + 0.05 * ( (1.-u) * vert4->x() + u * vert2->x() +
- (1.-v) * vert1->x() + v * vert3->x() -
- ( (1.-u)*(1.-v) * n1->x()
- + u * (1.-v) * n2->x()
- + u*v*n3->x()
- + (1.-u) * v * n4->x()) - vert->x());
- vert->y() = vert->y() + 0.05 * ( (1.-u) * vert4->y() + u * vert2->y() +
- (1.-v) * vert1->y() + v * vert3->y() -
- ( (1.-u)*(1.-v) * n1->y()
- + u * (1.-v) * n2->y()
- + u*v*n3->y()
- + (1.-u) * v * n4->y()) - vert->y());
+ Xold[i] = e[i]->x();
+ Yold[i] = e[i]->y();
+ Zold[i] = e[i]->z();
+ }
+
+ double minJ = 1.e22;
+ double maxJ = -1.e22;
+ getMinMaxJac (t1, minJ, maxJ);
+ getMinMaxJac (t2, minJ, maxJ);
+ int kopt = -1;
+ for (int k=0;k<NBST;k++){
+ double relax = (k+1)/(double)NBST;
+ for(unsigned int i = 0; i < e.size(); i++){
+ double v = (double)(i + 1) / (e.size() + 1);
+ double u = 1. - v;
+ MVertex *vert = (n2 < n4) ? e[i] : e[e.size() - i - 1];
+ MVertex *vert1 = (n1 < n2) ? e1[e1.size() - i - 1] : e1[i];
+ MVertex *vert3 = (n3 < n4) ? e3[i] : e3[e3.size() - i - 1];
+ MVertex *vert4 = (n4 < n1) ? e4[e4.size() - i - 1] : e4[i];
+ MVertex *vert2 = (n2 < n3) ? e2[i] : e2[e2.size() - i - 1];
+ double U1,V1,U2,V2,U3,V3,U4,V4,U,V,nU1,nV1,nU2,nV2,nU3,nV3,nU4,nV4;
+ parametricCoordinates(vert , gf,U,V);
+ parametricCoordinates(vert1, gf,U1,V1);
+ parametricCoordinates(vert2, gf,U2,V2);
+ parametricCoordinates(vert3, gf,U3,V3);
+ parametricCoordinates(vert4, gf,U4,V4);
+ parametricCoordinates(n1, gf,nU1,nV1);
+ parametricCoordinates(n2, gf,nU2,nV2);
+ parametricCoordinates(n3, gf,nU3,nV3);
+ parametricCoordinates(n4, gf,nU4,nV4);
+
+ Unew[k][i] = U + relax * ( (1.-u) * U4 + u * U2 +
+ (1.-v) * U1 + v * U3 -
+ ( (1.-u)*(1.-v) * nU1
+ + u * (1.-v) * nU2
+ + u*v*nU3
+ + (1.-u) * v * nU4) - U);
+ Vnew[k][i] = V + relax * ( (1.-u) * V4 + u * V2 +
+ (1.-v) * V1 + v * V3 -
+ ( (1.-u)*(1.-v) * nV1
+ + u * (1.-v) * nV2
+ + u*v*nV3
+ + (1.-u) * v * nV4) - V);
+ GPoint gp = gf->point(Unew[k][i],Vnew[k][i]);
+ vert->x() = gp.x();
+ vert->y() = gp.y();
+ vert->z() = gp.z();
+ }
+ double minJloc = 1.e22;
+ double maxJloc = -1.e22;
+ getMinMaxJac (t1, minJloc, maxJloc);
+ getMinMaxJac (t2, minJloc, maxJloc);
+
+ // printf("relax %d minJ %12.5E minjLoc %12.5E\n",k,minJ,minJloc);
+
+ if (minJloc > minJ){
+ kopt = k;
+ minJ = minJloc;
+ }
+ }
+ if (kopt == -1){
+ for(unsigned int i = 0; i < e.size(); i++){
+ e[i]->x() = Xold[i];
+ e[i]->y() = Yold[i];
+ e[i]->z() = Zold[i];
+ }
+ }
+ else{
+ success = true;
+ for(unsigned int i = 0; i < e.size(); i++){
+ MVertex *vert = (n2 < n4) ? e[i] : e[e.size() - i - 1];
+ vert->setParameter(0,Unew[kopt][i]);
+ vert->setParameter(1,Vnew[kopt][i]);
+ GPoint gp = gf->point(Unew[kopt][i],Vnew[kopt][i]);
+ vert->x() = gp.x();
+ vert->y() = gp.y();
+ vert->z() = gp.z();
+ }
}
}
}
}
+ return success;
}
void checkHighOrderTriangles(GModel *m)
{
+ double minJGlob = 1.e22;
+ double maxJGlob = -1.e22;
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
- bool twod = true;
+ double minJ = 1.e22;
+ double maxJ = -1.e22;
for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
+ double minJloc = 1.e22;
+ double maxJloc = -1.e22;
MTriangle *t = (*it)->triangles[i];
- if(t->getVertex(0)->z() != 0.0 ||
- t->getVertex(1)->z() != 0.0 ||
- t->getVertex(2)->z() != 0.0)
- twod = false;
- }
- if(twod){ // only perform the test for 2D/plane faces for now
- double minJ = 1.e22;
- double maxJ = -1.e22;
- for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
- double minJloc = 1.e22;
- double maxJloc = -1.e22;
- MTriangle *t = (*it)->triangles[i];
- if(t->getPolynomialOrder() > 1 && t->getPolynomialOrder() < 6){
- getMinMaxJac (t, minJloc, maxJloc);
- minJ = std::min(minJ, minJloc);
- maxJ = std::max(maxJ, maxJloc);
- if(minJloc * maxJloc < 0)
- Msg(WARNING, "Triangle %d (%d %d %d) has negative Jacobian", t->getNum(),
- t->getVertex(0)->getNum(), t->getVertex(1)->getNum(), t->getVertex(2)->getNum());
- }
- }
- if(minJ != 1.e22){
- if(minJ * maxJ < 0)
- Msg(GERROR, "Some triangles have negative Jacobians in surface %d", (*it)->tag());
- else
- Msg(INFO, "No negative Jacobians detected on model face %d: range = (%g,%g)",
- (*it)->tag(), minJ, maxJ);
+ if(t->getPolynomialOrder() > 1 && t->getPolynomialOrder() < 6){
+ getMinMaxJac (t, minJloc, maxJloc);
+ // printf("%p is %12.5E\n",t,minJloc);
+ minJ = std::min(minJ, minJloc);
+ maxJ = std::max(maxJ, maxJloc);
+// if(minJloc * maxJloc < 0)
+// Msg(WARNING, "Triangle %d (%d %d %d) has negative Jacobian (on gFace %d)", t->getNum(),
+// t->getVertex(0)->getNum(), t->getVertex(1)->getNum(), t->getVertex(2)->getNum(),(*it)->tag());
}
- }
+ }
+ minJGlob = std::min(minJGlob,minJ);
+ maxJGlob = std::max(maxJGlob,maxJ);
}
+ if (minJGlob >= 0)Msg(INFO, "Jacobian Range (%12.5E,%12.5E)",minJGlob, maxJGlob);
+ else Msg(WARNING, "Jacobian Range (%12.5E,%12.5E)",minJGlob, maxJGlob);
}
void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
@@ -866,11 +1186,28 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
if(CTX.mesh.smooth_internal_edges){
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
- for (int i = 0; i < 10; i++) smoothInternalEdges(*it, edgeVertices);
checkHighOrderTriangles(m);
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ Msg(INFO, "Smoothing internal Edges in Surface %d",(*it)->tag());
+ for (int i = 0; i < 10; i++) {
+ if (!smoothInternalEdges(*it, edgeVertices))break;
+ checkHighOrderTriangles(m);
+ }
+ }
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ for (int i=0;i<CTX.mesh.nb_smoothing;i++){
+ if(!optimizeHighOrderMesh(m,*it, edgeVertices))break;
+ checkHighOrderTriangles(m);
+ }
+ }
}
-
+
+ // MYFILE = fopen("jacs.pos","w");
+ // fprintf(MYFILE,"View \"\"{\n");
+ checkHighOrderTriangles(m);
+ // fprintf(MYFILE,"};\n");
+ // fclose(MYFILE);
+ MYFILE=0;
double t2 = Cpu();
Msg(INFO, "Mesh second order complete (%g s)", t2 - t1);
Msg(STATUS1, "Mesh");
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 3052065b43feb539e7fe1fabe1b196fdab016752..87c8e17b78016e5bd9b2295ce7a670405c7ce7e5 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGEdge.cpp,v 1.46 2007-11-11 19:53:57 remacle Exp $
+// $Id: meshGEdge.cpp,v 1.47 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -288,8 +288,6 @@ void meshGEdge::operator() (GEdge *ge)
if(ge->geomType() == GEntity::DiscreteCurve) return;
if(ge->geomType() == GEntity::BoundaryLayerCurve) return;
- // Send a messsage to the GMSH environment
- Msg(INFO, "Meshing curve %d", ge->tag());
deMeshGEdge dem;
dem(ge);
@@ -309,6 +307,10 @@ void meshGEdge::operator() (GEdge *ge)
// first compute the length of the curve by integrating one
double length = Integration(ge, t_begin, t_end, F_One, Points, 1.e-8);
ge->setLength(length);
+ // Send a messsage to the GMSH environment
+
+ if (length == 0.0)
+ Msg(GERROR,"Curve %d has a zero length\n",ge->tag());
List_Reset(Points);
@@ -337,6 +339,8 @@ void meshGEdge::operator() (GEdge *ge)
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.));
}
+ Msg(INFO, "Meshing curve %d (%s)", ge->tag(),ge->getTypeString().c_str());
+
// if the curve is periodic and if the begin vertex is identical to
// the end vertex and if this vertex has only one model curve
// adjacent to it, then the vertex is not connecting any other
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index d4e5cf8014973df9e4f02ea4a20ba08e37fc3a55..9358f946223df3c789e776d9745b3a886f92aa79 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.103 2008-01-09 15:25:48 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.104 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -35,11 +35,27 @@
#include "Numeric.h"
#include "BDS.h"
#include "qualityMeasures.h"
+#include "Field.h"
extern Context_T CTX;
static double SCALINGU=1,SCALINGV=1;
+bool noseam ( GFace *gf )
+{
+ std::list<GEdge*> edges = gf->edges();
+ std::list<GEdge*>::iterator it = edges.begin();
+ while (it != edges.end())
+ {
+ GEdge *ge = *it ;
+ bool seam = ge->isSeam(gf);
+ if (seam) return false;
+ ++it;
+ }
+ return true;
+}
+
+
void remeshUnrecoveredEdges ( std::set<EdgeToRecover> & edgesNotRecovered, std::list<GFace *> &facesToRemesh)
{
facesToRemesh.clear();
@@ -133,8 +149,8 @@ void remeshUnrecoveredEdges ( std::set<EdgeToRecover> & edgesNotRecovered, std::
bool AlgoDelaunay2D ( GFace *gf )
-{
- if ( /*gf->getNativeType() == GEntity::GmshModel &&*/ CTX.mesh.algo2d == ALGO_2D_DELAUNAY && gf->geomType() == GEntity::Plane)
+{
+ if ( noseam(gf) && /*gf->getNativeType() == GEntity::GmshModel &&*/ CTX.mesh.algo2d == ALGO_2D_DELAUNAY /*&& gf->geomType() == GEntity::Plane*/)
return true;
return false;
}
@@ -359,8 +375,19 @@ void OptimizeMesh(GFace *gf, BDS_Mesh &m, const int NIT)
}
}
}
- for (int KK=0;KK<4;KK++){
- // swap edges that provide a better configuration
+
+// for(int i = 0 ; i < NIT ; i++){
+// std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
+// while (itp != m.points.end())
+// {
+// optimize_vertex_position(gf,*itp,m.scalingU,m.scalingV);
+// ++itp;
+// }
+// }
+
+
+ for (int KK=0;KK<4;KK++){
+ // swap edges that provide a better configuration
int NN1 = m.edges.size();
int NN2 = 0;
std::list<BDS_Edge*>::iterator it = m.edges.begin();
@@ -398,7 +425,279 @@ void swapEdgePass ( GFace *gf, BDS_Mesh &m, int &nb_swap )
}
}
-void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
+int getTemplate (const BDS_Face *f, const std::map<BDS_Edge*,std::pair<BDS_Edge*,BDS_Edge*> >&splits,
+ BDS_Edge *edges[3][2]){
+ int k = 0;
+ std::map< BDS_Edge*,std::pair<BDS_Edge*,BDS_Edge* > > ::const_iterator it;
+ it = splits.find ( f-> e1 );
+ if (it == splits.end()){
+ edges[0][0] = f->e1;
+ edges[0][1] = 0;
+ }
+ else{
+ edges[0][0] = it->second.first;
+ edges[0][1] = it->second.second;
+ k+=1;
+ }
+ it = splits.find ( f->e2 );
+ if (it == splits.end()){
+ edges[1][0] = f->e2;
+ edges[1][1] = 0;
+ }
+ else{
+ edges[1][0] = it->second.first;
+ edges[1][1] = it->second.second;
+ k+=10;
+ }
+ it = splits.find ( f->e3 );
+ if (it == splits.end()){
+ edges[2][0] = f->e3;
+ edges[2][1] = 0;
+ }
+ else{
+ edges[2][0] = it->second.first;
+ edges[2][1] = it->second.second;
+ k+=100;
+ }
+ return k;
+}
+
+void Template_1 ( BDS_Mesh &m , BDS_Face *f, BDS_Edge *e11, BDS_Edge *e12, BDS_Edge *e2, BDS_Edge *e3)
+{
+ BDS_Point *mid = e11->commonvertex(e12);
+ BDS_Point *opposite = e2->commonvertex(e3);
+
+ if (!mid || !opposite){
+ printf("strange bazar in template 1 : edges %d %d , %d %d , %d %d and %d %d\n",
+ e11->p1->iD,e11->p2->iD,
+ e12->p1->iD,e12->p2->iD,
+ e2->p1->iD,e2->p2->iD,
+ e3->p1->iD,e3->p2->iD);
+ }
+
+ BDS_Edge *emid = new BDS_Edge (mid,opposite);
+
+ if (!e11->commonvertex(e3)){
+ BDS_Edge *temp = e3;
+ e3 = e2;
+ e2 = temp;
+ }
+
+ BDS_Face* t1 = new BDS_Face(e11,emid,e3);
+ BDS_Face* t2 = new BDS_Face(e12,e2,emid);
+ t1->g = f->g;
+ t2->g = t1->g;
+ emid->g = t1->g;
+ m.triangles.push_back(t1);
+ m.triangles.push_back(t2);
+ m.edges.push_back(emid);
+ m.del_face(f);
+ e11->p1->config_modified = e11->p2->config_modified = true;
+ e12->p1->config_modified = e12->p2->config_modified = true;
+ e2->p1->config_modified = e2->p2->config_modified = true;
+ e3->p1->config_modified = e3->p2->config_modified = true;
+}
+
+void Template_2 ( BDS_Mesh &m , BDS_Face *f, BDS_Edge *e11, BDS_Edge *e12, BDS_Edge *e21, BDS_Edge *e22, BDS_Edge *e3)
+{
+ BDS_Point *mid1 = e11->commonvertex(e12);
+ BDS_Point *mid2 = e21->commonvertex(e22);
+
+ BDS_Edge *emid1 = new BDS_Edge (mid1,mid2);
+
+ if (!e11->commonvertex(e3)){
+ BDS_Edge *temp = e11;
+ e11 = e12;
+ e12 = temp;
+ }
+ if (!e22->commonvertex(e3)){
+ BDS_Edge *temp = e22;
+ e22 = e21;
+ e21 = temp;
+ }
+ BDS_Point *opposite1 = e3->commonvertex(e22);
+ BDS_Point *opposite2 = e3->commonvertex(e11);
+
+ // if (!e12->commonvertex(e21))throw;
+
+ // build the best possible template to avoid subsequent swap
+ // first config, use and edge mid1->opposite1
+ double config1Q = std::min(qmTriangle(opposite2,mid1,opposite1,QMTRI_RHO),qmTriangle(mid1,mid2,opposite1,QMTRI_RHO));
+ // second config, use and edge mid2->opposite2
+ double config2Q = std::min(qmTriangle(opposite2,mid2,opposite1,QMTRI_RHO),qmTriangle(mid1,mid2,opposite2,QMTRI_RHO));
+
+ // if the first one is the best
+ BDS_Face *t1,*t2,*t3;
+ BDS_Edge *emid2;
+ t1 = new BDS_Face(e12,e21,emid1);
+ if (config1Q > config2Q) {
+ emid2 = new BDS_Edge (mid1,opposite1);
+ t2 = new BDS_Face(e11,emid2,e3);
+ t3 = new BDS_Face(emid2,emid1,e22);
+ }
+ else{
+ emid2 = new BDS_Edge (mid2,opposite2);
+ t2 = new BDS_Face(e11,emid1,emid2);
+ t3 = new BDS_Face(emid2,e22,e3);
+ }
+ t1->g = f->g;
+ t2->g = t1->g;
+ t3->g = t1->g;
+ emid1->g = t1->g;
+ emid2->g = t1->g;
+ m.triangles.push_back(t1);
+ m.triangles.push_back(t2);
+ m.triangles.push_back(t3);
+ m.edges.push_back(emid1);
+ m.edges.push_back(emid2);
+ m.del_face(f);
+ e11->p1->config_modified = e11->p2->config_modified = true;
+ e12->p1->config_modified = e12->p2->config_modified = true;
+ e21->p1->config_modified = e21->p2->config_modified = true;
+ e22->p1->config_modified = e22->p2->config_modified = true;
+ e3->p1->config_modified = e3->p2->config_modified = true;
+}
+void Template_3 ( BDS_Mesh &m , BDS_Face *f, BDS_Edge *e11, BDS_Edge *e12, BDS_Edge *e21, BDS_Edge *e22, BDS_Edge *e31,BDS_Edge *e32)
+{
+ BDS_Point *mid1 = e11->commonvertex(e12);
+ BDS_Point *mid2 = e21->commonvertex(e22);
+ BDS_Point *mid3 = e31->commonvertex(e32);
+
+ BDS_Edge *emid12 = new BDS_Edge (mid1,mid2);
+ BDS_Edge *emid13 = new BDS_Edge (mid1,mid3);
+ BDS_Edge *emid23 = new BDS_Edge (mid2,mid3);
+
+ if (!e11->commonvertex(e31) && !e11->commonvertex(e32)){
+ BDS_Edge *temp = e11;
+ e11 = e12;
+ e12 = temp;
+ }
+ if (!e11->commonvertex(e31)){
+ BDS_Edge *temp = e31;
+ e31 = e32;
+ e32 = temp;
+ }
+ if (!e32->commonvertex(e22)){
+ BDS_Edge *temp = e21;
+ e21 = e22;
+ e22 = temp;
+ }
+
+// if (!e11->commonvertex(e31))throw;
+// if (!e32->commonvertex(e22))throw;
+// if (!e12->commonvertex(e21))throw;
+
+ BDS_Face *t1 = new BDS_Face(e11,emid13,e31);
+ BDS_Face *t2 = new BDS_Face(e12,e21,emid12);
+ BDS_Face *t3 = new BDS_Face(emid23,e22,e32);
+ BDS_Face *t4 = new BDS_Face(emid12,emid23,emid13);
+
+ t1->g = f->g;
+ t2->g = t1->g;
+ t3->g = t1->g;
+ t4->g = t1->g;
+ emid12->g = t1->g;
+ emid13->g = t1->g;
+ emid23->g = t1->g;
+ m.triangles.push_back(t1);
+ m.triangles.push_back(t2);
+ m.triangles.push_back(t3);
+ m.triangles.push_back(t4);
+ m.edges.push_back(emid12);
+ m.edges.push_back(emid13);
+ m.edges.push_back(emid23);
+ m.del_face(f);
+ e11->p1->config_modified = e11->p2->config_modified = true;
+ e12->p1->config_modified = e12->p2->config_modified = true;
+ e21->p1->config_modified = e21->p2->config_modified = true;
+ e22->p1->config_modified = e22->p2->config_modified = true;
+ e31->p1->config_modified = e31->p2->config_modified = true;
+ e32->p1->config_modified = e32->p2->config_modified = true;
+}
+
+
+void splitEdgePass_templateRefine ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
+{
+ std::map<BDS_Edge*,std::pair<BDS_Edge*,BDS_Edge*> >splits;
+ // build a list of edges to split with their new vertices
+ int NN1 = m.edges.size();
+ int NN2 = 0;
+ std::list<BDS_Edge*>::iterator it = m.edges.begin();
+ while (1)
+ {
+ if (NN2++ >= NN1)break;
+ if (!(*it)->deleted && (*it)->numfaces() == 2)
+ {
+ double lone = NewGetLc ( *it,gf);
+ if (lone > MAXE_){
+ const double coord = 0.5;
+ BDS_Point *mid ;
+ mid = m.add_point(++m.MAXPOINTNUMBER,
+ coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
+ coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v,gf);
+
+ mid->lcBGM() = BGM_MeshSize(gf,
+ (coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u)*m.scalingU,
+ (coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
+ mid->X,mid->Y,mid->Z);
+ mid->lc() = 0.5 * ( (*it)->p1->lc() + (*it)->p2->lc() );
+ BDS_Edge *e1 = new BDS_Edge ((*it)->p1,mid);
+ BDS_Edge *e2 = new BDS_Edge (mid, (*it)->p2);
+ e1->g = e2->g = (*it)->g;
+ m.del_edge((*it));
+ m.edges.push_back(e1);
+ m.edges.push_back(e2);
+
+
+ splits[*it] = std::make_pair<BDS_Edge*, BDS_Edge*> (e1,e2);
+ nb_split++;
+ }
+ }
+ ++it;
+ }
+
+ std::list<BDS_Face*>::iterator itt = m.triangles.begin();
+ // build a list of edges to split with their new vertices
+ while (itt != m.triangles.end())
+ {
+ if (!(*itt)->deleted)
+ {
+ BDS_Edge *edges[3][2];
+ int K = getTemplate ((*itt),splits,edges);
+ switch(K){
+ case 0: // no edge is split
+ break;
+ case 1: // first edge is split
+ Template_1 ( m , *itt, edges[0][0],edges[0][1],edges[1][0],edges[2][0]);
+ break;
+ case 10: // second edge is split
+ Template_1 ( m , *itt, edges[1][0],edges[1][1],edges[2][0],edges[0][0]);
+ break;
+ case 100: // third edge is split
+ Template_1 ( m , *itt, edges[2][0],edges[2][1],edges[0][0],edges[1][0]);
+ break;
+ case 11: // fisrt and second
+ Template_2 ( m , *itt, edges[0][0],edges[0][1],edges[1][0],edges[1][1],edges[2][0]);
+ break;
+ case 101: // fisrt and third
+ Template_2 ( m , *itt, edges[2][0],edges[2][1],edges[0][0],edges[0][1],edges[1][0]);
+ break;
+ case 110: // second and third
+ Template_2 ( m , *itt, edges[1][0],edges[1][1],edges[2][0],edges[2][1],edges[0][0]);
+ break;
+ case 111: // all splitted
+ Template_3 ( m , *itt, edges[0][0],edges[0][1],edges[1][0],edges[1][1],edges[2][0],edges[2][1]);
+ break;
+ default :
+ printf("strange template %d\n",K);
+ throw;
+ }
+ }
+ ++itt;
+ }
+}
+
+void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split, AttractorField_1DMesh *attr)
{
int NN1 = m.edges.size();
int NN2 = 0;
@@ -428,7 +727,21 @@ void splitEdgePass ( GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
(coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v)*m.scalingV,
mid->X,mid->Y,mid->Z);
//mid->lc() = 2./ ( 1./(*it)->p1->lc() + 1./(*it)->p2->lc() );
- mid->lc() = 0.5 * ( (*it)->p1->lc() + (*it)->p2->lc() );
+ if (!attr)
+ mid->lc() = 0.5 * ( (*it)->p1->lc() + (*it)->p2->lc() );
+ else{
+ double lcmin,lcpt, dist;
+ double FACT1 = 3, FACT2=25;
+ attr->eval(mid->X,mid->Y,mid->Z,lcmin,lcpt, dist);
+
+ if (dist < FACT1* lcmin) mid->lc() = lcmin;
+ else if (dist < FACT2 * lcmin) {
+ double r = (dist - FACT1*lcmin) / ((FACT2-FACT1)*lcmin);
+ mid->lc() = 1./(1./lcmin * r + 1./lcpt*(1-r));
+ }
+ else mid->lc() = lcpt;
+ }
+
// printf("%g %g\n",mid->lc(),mid->lcBGM());
if(!m.split_edge ( *it, mid )) m.del_point(mid);
@@ -598,51 +911,45 @@ void smoothVertexPass ( GFace *gf, BDS_Mesh &m, int &nb_smooth)
void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
{
int IT =0;
+ AttractorField_1DMesh *attr = 0;
+// if (CTX.mesh.lc_from_curvature){
+// // parameters are not important
+// attr = new AttractorField_1DMesh(gf,0.1,0,1);
+// attr->buildFastSearchStructures();
+// }
// printf("lc (1,1) = %g\n",Attractor::lc(1,1,0));
int MAXNP = m.MAXPOINTNUMBER;
- // if (NIT > 0)
+ if (NIT > 0)
{
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
while (itp != m.points.end())
{
- if (NIT > 0)
- {
- std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
- std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
- double L=0;
- int ne = 0;
- while(it!=ite){
- double l = (*it)->length();
- if ((*it)->g && (*it)->g->classif_degree == 1){
- L=ne?std::max(L,l):l;
- // L=ne?std::min(L,l):l;
- // L+=l;
- ne++;
- }
- ++it;
- }
- if (!ne) L = 1.e22;
- // else L/=ne;
- if(!CTX.mesh.constrained_bgmesh)
- (*itp)->lc() = L;
- (*itp)->lcBGM() = L;
+ std::list<BDS_Edge*>::iterator it = (*itp)->edges.begin();
+ std::list<BDS_Edge*>::iterator ite = (*itp)->edges.end();
+ double L=0;
+ int ne = 0;
+ while(it!=ite){
+ double l = (*it)->length();
+ if ((*it)->g && (*it)->g->classif_degree == 1){
+ L=ne?std::max(L,l):l;
+ // L=ne?std::min(L,l):l;
+ // L+=l;
+ ne++;
}
-// else if (CTX.mesh.lc_from_curvature)
-// {
-// double Crv = gf->curvature(SPoint2((*itp)->u, (*itp)->v));
-// double lc = Crv > 0 ? 2*M_PI / Crv / CTX.mesh.min_circ_points : 1.e22;
-// lc *= CTX.mesh.lc_factor;
-// printf("%d = %d %12.5e %12.5E %12.5E\n",gf->tag(),CTX.mesh.min_circ_points,Crv,lc,(*itp)->lcBGM());
-// (*itp)->lc() = std::min(lc,(*itp)->lc());
-// }
+ ++it;
+ }
+ if (!ne) L = 1.e22;
+ // else L/=ne;
+ if(!CTX.mesh.constrained_bgmesh)
+ (*itp)->lc() = L;
+ (*itp)->lcBGM() = L;
++itp;
}
}
-
double OLDminL=1.E22,OLDmaxL=0;
double t_spl=0, t_sw=0,t_col=0,t_sm=0;
@@ -684,11 +991,12 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
OLDminL = minL;OLDmaxL = maxL;
if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT)))break;
-
- double maxE = MAXE_;//std::max(MAXE_,maxL / 5);
- double minE = MINE_;//std::min(MINE_,minL * 1.2);
+ double maxE = MAXE_;
+ double minE = MINE_;
clock_t t1 = clock();
- splitEdgePass ( gf, m, maxE, nb_split);
+ // splitEdgePass_templateRefine ( gf, m, maxE, nb_split);
+ splitEdgePass ( gf, m, maxE, nb_split,attr);
+ //splitEdgePass_templateRefine ( gf, m, maxE, nb_split);
//saturateEdgePass ( gf, m, maxE, nb_split);
clock_t t2 = clock();
swapEdgePass ( gf, m, nb_swap);
@@ -721,6 +1029,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
if (nb_split==0 && nb_collaps == 0)break;
// if (mesh_quality < old_mesh_quality && worst > 0.1 && maxL < 1.45) break;
+ // return;
+
}
double t_total = t_spl + t_sw + t_col + t_sm;
@@ -734,6 +1044,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
Msg(DEBUG1," ---------------------------------------");
Msg(DEBUG1," CPU TOTAL %12.5E sec ",t_total);
Msg(DEBUG1," ---------------------------------------");
+
+ if (attr)delete attr;
}
@@ -1196,9 +1508,6 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
}
m->cleanup();
-
-
-
{
std::list<BDS_Edge*>::iterator ite = m->edges.begin();
while (ite != m->edges.end())
@@ -1225,7 +1534,6 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
m->del_point(m->find_point(-3));
m->del_point(m->find_point(-4));
-
if (debug)
{
char name[245];
@@ -1238,10 +1546,11 @@ bool gmsh2DMeshGenerator ( GFace *gf , int RECUR_ITER, bool debug = true)
// start mesh generation
if (!AlgoDelaunay2D ( gf ))
{
+ // RefineMesh (gf,*m, 1);
RefineMesh (gf,*m, CTX.mesh.refine_steps);
OptimizeMesh(gf, *m, 2);
RefineMesh (gf,*m, -CTX.mesh.refine_steps);
- OptimizeMesh(gf, *m, 2);
+ OptimizeMesh(gf, *m, -2);
if (gf->meshAttributes.recombine)
{
@@ -1341,20 +1650,6 @@ inline double dist2 (const SPoint2 &p1,const SPoint2 &p2)
-bool noseam ( GFace *gf )
-{
- std::list<GEdge*> edges = gf->edges();
- std::list<GEdge*>::iterator it = edges.begin();
- while (it != edges.end())
- {
- GEdge *ge = *it ;
- bool seam = ge->isSeam(gf);
- if (seam) return false;
- ++it;
- }
- return true;
-}
-
bool buildConsecutiveListOfVertices ( GFace *gf,
GEdgeLoop &gel ,
std::vector<BDS_Point*> &result,
@@ -1698,16 +1993,14 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
BDS_GeomEntity CLASS_F (1,2);
BDS_GeomEntity CLASS_E (1,1);
-
- if (debug)
- {
- char name[245];
- sprintf(name,"surface%d-initial-real.pos",gf->tag());
- outputScalarField(m->triangles, name,0);
- sprintf(name,"surface%d-initial-param.pos",gf->tag());
- outputScalarField(m->triangles, name,1);
- }
+ if (debug){
+ char name[245];
+ sprintf(name,"surface%d-initial-real.pos",gf->tag());
+ outputScalarField(m->triangles, name,0);
+ sprintf(name,"surface%d-initial-param.pos",gf->tag());
+ outputScalarField(m->triangles, name,1);
+ }
// build a list of edges to recover
@@ -1722,22 +2015,19 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
// }
- for (unsigned int i=0;i<edgeLoops_BDS.size();i++)
- {
+ for (unsigned int i=0;i<edgeLoops_BDS.size();i++){
std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
- for (unsigned int j=0;j<edgeLoop_BDS.size();j++)
- {
- BDS_Edge * e = m->recover_edge ( edgeLoop_BDS[j]->iD,edgeLoop_BDS[(j+1)%edgeLoop_BDS.size()]->iD);
- if (!e)
- {
- Msg(GERROR,"impossible to recover the edge %d %d",edgeLoop_BDS[j]->iD,edgeLoop_BDS[(j+1)%edgeLoop_BDS.size()]->iD);
- gf->meshStatistics.status = GFace::FAILED;
- SCALINGU = SCALINGV = 1;
- return false;
- }
- else e->g = &CLASS_E;
+ for (unsigned int j=0;j<edgeLoop_BDS.size();j++){
+ BDS_Edge * e = m->recover_edge ( edgeLoop_BDS[j]->iD,edgeLoop_BDS[(j+1)%edgeLoop_BDS.size()]->iD);
+ if (!e){
+ Msg(GERROR,"impossible to recover the edge %d %d",edgeLoop_BDS[j]->iD,edgeLoop_BDS[(j+1)%edgeLoop_BDS.size()]->iD);
+ gf->meshStatistics.status = GFace::FAILED;
+ SCALINGU = SCALINGV = 1;
+ return false;
}
- }
+ else e->g = &CLASS_E;
+ }
+ }
// Msg(INFO,"Boundary Edges recovered for surface %d",gf->tag());
// Look for an edge that is on the boundary for which one of the
@@ -1747,40 +2037,34 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
// recursive algorithm
{
std::list<BDS_Edge*>::iterator ite = m->edges.begin();
- while (ite != m->edges.end())
- {
- BDS_Edge *e = *ite;
- if ( e->g && e->numfaces () == 2)
- {
- BDS_Point *oface[2];
- e->oppositeof(oface);
- if (oface[0]->iD < 0)
- {
- recur_tag ( e->faces(1) , &CLASS_F);
- break;
- }
- else if (oface[1]->iD < 0)
- {
- recur_tag ( e->faces(0) , &CLASS_F);
- break;
- }
- }
- ++ite;
+ while (ite != m->edges.end()){
+ BDS_Edge *e = *ite;
+ if ( e->g && e->numfaces () == 2){
+ BDS_Point *oface[2];
+ e->oppositeof(oface);
+ if (oface[0]->iD < 0){
+ recur_tag ( e->faces(1) , &CLASS_F);
+ break;
+ }
+ else if (oface[1]->iD < 0){
+ recur_tag ( e->faces(0) , &CLASS_F);
+ break;
+ }
}
+ ++ite;
+ }
}
// delete useless stuff
{
std::list<BDS_Face*>::iterator itt = m->triangles.begin();
- while (itt != m->triangles.end())
- {
- BDS_Face *t = *itt;
- if (!t->g)
- {
- m->del_face (t);
- }
- ++itt;
+ while (itt != m->triangles.end()){
+ BDS_Face *t = *itt;
+ if (!t->g){
+ m->del_face (t);
}
+ ++itt;
+ }
}
m->cleanup();
@@ -1826,7 +2110,7 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
RefineMesh (gf,*m,CTX.mesh.refine_steps);
OptimizeMesh(gf, *m, 2);
RefineMesh (gf,*m,-CTX.mesh.refine_steps);
- OptimizeMesh(gf, *m, 2);
+ OptimizeMesh(gf, *m, -2);
if (gf->meshAttributes.recombine)
{
@@ -1965,10 +2249,11 @@ void meshGFace::operator() (GFace *gf)
computeEdgeLoops(gf, points, indices);
// temp fix until we create MEdgeLoops in gmshFace
- // if (gf->tag() == 2)
+ if (1 || gf->tag() == 6)
{
Msg(DEBUG1, "Generating the mesh");
if(noseam (gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty()){
+ // gmsh2DMeshGenerator(gf,0, true);
gmsh2DMeshGenerator(gf,0, false);
}
else{
@@ -1976,9 +2261,9 @@ void meshGFace::operator() (GFace *gf)
Msg(GERROR, "Impossible to mesh face %d", gf->tag());
}
}
-// else
-// gf->meshStatistics.status = GFace::DONE;
-
+ else
+ gf->meshStatistics.status = GFace::DONE;
+
Msg(DEBUG1, "type %d %d triangles generated, %d internal vertices",
gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size());
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 99b1cadd6c689dd4b702a7e947a20b775773a637..42cd02777dd5b7b9ab4e752f57b057608179591c 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.6 2007-10-10 08:49:34 remacle Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.7 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -22,6 +22,7 @@
#include "BDS.h"
#include "BackgroundMesh.h"
#include "meshGFaceDelaunayInsertion.h"
+#include "meshGFaceOptimize.h"
#include "meshGFace.h"
#include "GFace.h"
#include "Numeric.h"
@@ -32,41 +33,100 @@
int III = 1;
-// computes the center of the circum circle in the tangent plane
-// the metric given by a b d is supposed to be constant
-void MTri3::Center_Circum_Aniso(double a, double b, double d, double &x, double &y, double &r) const
+bool circumCenterMetricInTriangle ( MTriangle *base,
+ const double *metric,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs)
{
- throw;
+ double R,x[2],uv[2];
+ circumCenterMetric(base,metric,Us,Vs,x,R);
+ return invMapUV(base,x,Us,Vs,uv,1.e-8);
}
-
-// find a common basis for 2 metrics in 2D
-void simultaneousMetricReduction( const gmsh2dMetric &M1, const gmsh2dMetric &M2,
- double & l1,double & l2, double V[2][2])
+void circumCenterMetric ( MTriangle *base,
+ const double *metric,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs,
+ double *x, double &Radius2)
{
- throw;
-}
+ // d = (u2-u1) M (u2-u1) = u2 M u2 + u1 M u1 - 2 u2 M u1
+ double sys[2][2];
+ double rhs[2];
-gmsh2dMetric metricIntersection(const gmsh2dMetric &Ma,const gmsh2dMetric &Mb)
-{
- throw;
+ double pa[2] = {Us[base->getVertex(0)->getNum()],
+ Vs[base->getVertex(0)->getNum()]};
+ double pb[2] = {Us[base->getVertex(1)->getNum()],
+ Vs[base->getVertex(1)->getNum()]};
+ double pc[2] = {Us[base->getVertex(2)->getNum()],
+ Vs[base->getVertex(2)->getNum()]};
+
+ const double a = metric[0];
+ const double b = metric[1];
+ const double d = metric[2];
+
+ sys[0][0] = 2. * a * (pa[0] - pb[0]) + 2. * b * (pa[1] - pb[1]);
+ sys[0][1] = 2. * d * (pa[1] - pb[1]) + 2. * b * (pa[0] - pb[0]);
+ sys[1][0] = 2. * a * (pa[0] - pc[0]) + 2. * b * (pa[1] - pc[1]);
+ sys[1][1] = 2. * d * (pa[1] - pc[1]) + 2. * b * (pa[0] - pc[0]);
+
+ rhs[0] =
+ a * (pa[0] * pa[0] - pb[0] * pb[0]) +
+ d * (pa[1] * pa[1] - pb[1] * pb[1]) +
+ 2. * b * (pa[0] * pa[1] - pb[0] * pb[1]);
+ rhs[1] =
+ a * (pa[0] * pa[0] - pc[0] * pc[0]) +
+ d * (pa[1] * pa[1] - pc[1] * pc[1]) +
+ 2. * b * (pa[0] * pa[1] - pc[0] * pc[1]);
+
+ sys2x2(sys, rhs, x);
+
+ Radius2 =
+ (x[0] - pa[0]) * (x[0] - pa[0]) * a +
+ (x[1] - pa[1]) * (x[1] - pa[1]) * d +
+ 2. * (x[0] - pa[0]) * (x[1] - pa[1]) * b;
}
-gmsh2dMetric metricInterpolationTriangle(const gmsh2dMetric &M1,
- const gmsh2dMetric &M2,
- const gmsh2dMetric &M3,
- const double u,
- const double v)
+void buildMetric ( GFace *gf , double *uv, double *metric){
+ Pair<SVector3,SVector3> der = gf->firstDer(SPoint2(uv[0],uv[1]));
+ metric[0] = dot(der.first(),der.first());
+ metric[1] = dot(der.second(),der.first());
+ metric[2] = dot(der.second(),der.second());
+}
+
+
+int inCircumCircleAniso ( GFace *gf,
+ MTriangle *base,
+ const double *uv ,
+ const double *metricb,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs)
{
- throw;
-}
-gmsh2dMetric :: gmsh2dMetric ( double lc )
- : a11 ( 1./(lc*lc) ) ,a21 ( 0.0 ) ,a22 ( 1./(lc*lc) )
-{
- throw;
-}
+ double x[2],Radius2,metric[3];
+ double pa[2] = {(Us[base->getVertex(0)->getNum()]+Us[base->getVertex(1)->getNum()]+Us[base->getVertex(2)->getNum()])/3.,
+ (Vs[base->getVertex(0)->getNum()]+Vs[base->getVertex(1)->getNum()]+Vs[base->getVertex(2)->getNum()])/3.};
+
+ buildMetric (gf,pa,metric);
+ circumCenterMetric ( base, metric, Us,Vs,x,Radius2);
+// buildMetric (gf,pa,metric);
+// circumCenterMetric ( base, metric, Us,Vs,x,Radius2);
+// buildMetric (gf,pa,metric);
+// circumCenterMetric ( base, metric, Us,Vs,x,Radius2);
+// buildMetric (gf,pa,metric);
+// circumCenterMetric ( base, metric, Us,Vs,x,Radius2);
+
+ const double a = metric[0];
+ const double b = metric[1];
+ const double d = metric[2];
+
+
+ double d2 = (x[0] - uv[0]) * (x[0] - uv[0]) * a
+ + (x[1] - uv[1]) * (x[1] - uv[1]) * d
+ + 2. * (x[0] - uv[0]) * (x[1] - uv[1]) * b;
+
+ return d2 < Radius2;
+}
MTri3::MTri3 ( MTriangle * t, double lc) : deleted(false), base (t)
{
@@ -159,6 +219,10 @@ void connectTriangles ( std::list<MTri3*> & l)
{
connectTris(l.begin(),l.end());
}
+void connectTriangles ( std::vector<MTri3*> & l)
+{
+ connectTris(l.begin(),l.end());
+}
void recurFindCavity (std::list<edgeXface> & shell,
std::list<MTri3*> & cavity,
@@ -189,6 +253,38 @@ void recurFindCavity (std::list<edgeXface> & shell,
}
}
+void recurFindCavityAniso (GFace *gf,
+ std::list<edgeXface> & shell,
+ std::list<MTri3*> & cavity,
+ double *metric,
+ double *param,
+ MTri3 *t,
+ std::vector<double> & Us,
+ std::vector<double> & Vs)
+{
+ t->setDeleted(true);
+ // the cavity that has to be removed
+ // because it violates delaunay criterion
+ cavity.push_back(t);
+
+ for (int i=0;i<3;i++)
+ {
+ MTri3 *neigh = t->getNeigh(i) ;
+ if (!neigh)
+ shell.push_back ( edgeXface ( t, i ) );
+ else if (!neigh->isDeleted())
+ {
+ int circ = inCircumCircleAniso (gf, neigh->tri(), param, metric, Us, Vs);
+ if (circ)
+ recurFindCavityAniso ( gf, shell, cavity,metric, param, neigh,Us,Vs);
+ else
+ shell.push_back ( edgeXface ( t, i ) );
+ }
+ }
+}
+
+
+
bool circUV ( MTriangle *t ,
std::vector<double> & Us,
std::vector<double> & Vs , double *res, GFace *gf)
@@ -214,8 +310,8 @@ bool circUV ( MTriangle *t ,
bool invMapUV ( MTriangle *t ,
double *p,
- std::vector<double> & Us,
- std::vector<double> & Vs , double *uv, double tol)
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs , double *uv, double tol)
{
double mat[2][2];
double b[2];
@@ -264,23 +360,29 @@ double getSurfUV ( MTriangle *t ,
}
-bool insertVertex (MVertex *v ,
- double *param ,
+bool insertVertex (GFace *gf,
+ MVertex *v ,
+ double *param ,
MTri3 *t ,
std::set<MTri3*,compareTri3Ptr> &allTets,
std::vector<double> & vSizes,
std::vector<double> & vSizesBGM,
std::vector<double> & Us,
- std::vector<double> & Vs)
+ std::vector<double> & Vs,
+ double *metric = 0)
{
std::list<edgeXface> shell;
std::list<MTri3*> cavity;
std::list<MTri3*> new_cavity;
- double p[3] = {v->x(),v->y(),v->z()};
-
- recurFindCavity ( shell, cavity, p , param, t, Us, Vs);
-
+ if (!metric){
+ double p[3] = {v->x(),v->y(),v->z()};
+ recurFindCavity ( shell, cavity, p , param, t, Us, Vs);
+ }
+ else{
+ recurFindCavityAniso ( gf, shell, cavity, metric , param, t, Us, Vs);
+ }
+
// check that volume is conserved
double newVolume = 0;
double oldVolume = 0;
@@ -361,7 +463,7 @@ bool insertVertex (MVertex *v ,
// getchar();
- if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume )
+ if (fabs(oldVolume - newVolume) < 1.e-12 * oldVolume )
{
connectTris ( new_cavity.begin(),new_cavity.end() );
allTets.insert(newTris,newTris+shell.size());
@@ -482,6 +584,13 @@ void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
gf->triangles.clear();
connectTris ( AllTris.begin(), AllTris.end() );
+
+ // _printTris ("before.pos", AllTris, Us,Vs);
+ // this should be MUCH faster !
+ for (int i=0;i<1200;i++)
+ if(!edgeSwapPass(gf,AllTris,SWCR_DEL,Us,Vs,vSizes,vSizesBGM))break;
+
+ // _printTris ("after2.pos", AllTris, Us,Vs);
Msg(DEBUG,"All %d tris were connected",AllTris.size());
@@ -493,90 +602,111 @@ void insertVerticesInFace (GFace *gf, BDS_Mesh *bds)
{
MTri3 *worst = *AllTris.begin();
- if (worst->isDeleted())
- {
- delete worst->tri();
- delete worst;
- AllTris.erase(AllTris.begin());
- // Msg(INFO,"Worst tet is deleted");
+ if (worst->isDeleted()){
+ delete worst->tri();
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ // Msg(INFO,"Worst tet is deleted");
+ }
+ else{
+ if(ITER++%5000 ==0)
+ Msg(DEBUG,"%7d points created -- Worst tri radius is %8.3f",vSizes.size(),worst->getRadius());
+ double center[2],uv[2],metric[3],r2;
+ if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
+
+ circUV(worst->tri(),Us,Vs,center,gf);
+ MTriangle *base = worst->tri();
+ double pa[2] = {(Us[base->getVertex(0)->getNum()]+Us[base->getVertex(1)->getNum()]+Us[base->getVertex(2)->getNum()])/3.,
+ (Vs[base->getVertex(0)->getNum()]+Vs[base->getVertex(1)->getNum()]+Vs[base->getVertex(2)->getNum()])/3.};
+ buildMetric ( gf , pa, metric);
+ circumCenterMetric ( worst->tri(),metric,Us,Vs,center,r2);
+ // buildMetric ( gf , center, metric);
+ // circumCenterMetric ( worst->tri(),metric,Us,Vs,center,r2);
+ // buildMetric ( gf , center, metric);
+ // circumCenterMetric ( worst->tri(),metric,Us,Vs,center,r2);
+
+ // printf("%g %g %g\n",metric[0],metric[1],metric[2]);
+
+ bool inside = invMapUV(worst->tri(),center,Us,Vs,uv,1.e-8);
+ // if (!inside)
+ // circUV(worst->tri(),Us,Vs,center,gf);
+
+ // inside = invMapUV(worst->tri(),center,Us,Vs,uv,1.e-8);
+
+ if (inside) {
+ // we use here local coordinates as real coordinates
+ // x,y and z will be computed hereafter
+ // Msg(INFO,"Point is inside");
+ GPoint p = gf->point (center[0],center[1]);
+ MVertex *v = new MFaceVertex (p.x(),p.y(),p.z(),gf,center[0],center[1]);
+ v->setNum(NUM++);
+ double lc1 = ((1.-uv[0]-uv[1]) * vSizes [worst->tri()->getVertex(0)->getNum()] +
+ uv[0] * vSizes [worst->tri()->getVertex(1)->getNum()] +
+ uv[1] * vSizes [worst->tri()->getVertex(2)->getNum()] );
+ // double eigMetricSurface = gf->getMetricEigenvalue(SPoint2(center[0],center[1]));
+ double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
+ // printf("lc1 %12.5E lc %12.5E\n",lc1,lc);
+
+ vSizesBGM.push_back( lc );
+ vSizes.push_back ( lc1);
+ Us.push_back( center[0] );
+ Vs.push_back( center[1] );
+
+ if (!insertVertex ( gf, v , center, worst, AllTris,vSizes,vSizesBGM,Us,Vs,metric)) {
+ Msg(DEBUG,"2D Delaunay : a cavity is not star shaped");
+ AllTris.erase(AllTris.begin());
+ worst->forceRadius(-1);
+ AllTris.insert(worst);
+ delete v;
+ }
+ else
+ gf->mesh_vertices.push_back(v);
}
- else
- {
- if(ITER++%5000 ==0)
- Msg(DEBUG,"%7d points created -- Worst tri radius is %8.3f",vSizes.size(),worst->getRadius());
- double center[2],uv[2];
- if (worst->getRadius() < 0.5 * sqrt(2.0)) break;
- circUV(worst->tri(),Us,Vs,center,gf);
- bool inside = invMapUV(worst->tri(),center,Us,Vs,uv,1.e-8);
-
- if (inside)
- {
-
-// char name[245];
-// sprintf(name,"param%d.pos",ITER);
-// _printTris (name, AllTris, Us,Vs);
-
-
- // we use here local coordinates as real coordinates
- // x,y and z will be computed hereafter
- // Msg(INFO,"Point is inside");
- GPoint p = gf->point (center[0],center[1]);
- MVertex *v = new MFaceVertex (p.x(),p.y(),p.z(),gf,center[0],center[1]);
- v->setNum(NUM++);
- double lc1 = ((1.-uv[0]-uv[1]) * vSizes [worst->tri()->getVertex(0)->getNum()] +
- uv[0] * vSizes [worst->tri()->getVertex(1)->getNum()] +
- uv[1] * vSizes [worst->tri()->getVertex(2)->getNum()] );
- // double eigMetricSurface = gf->getMetricEigenvalue(SPoint2(center[0],center[1]));
- double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
- // printf("lc1 %12.5E lc %12.5E\n",lc1,lc);
-
- vSizesBGM.push_back( lc );
- vSizes.push_back ( lc1);
- Us.push_back( center[0] );
- Vs.push_back( center[1] );
-
- if (!insertVertex ( v , center, worst, AllTris,vSizes,vSizesBGM,Us,Vs))
- {
- Msg(DEBUG,"2D Delaunay : a cavity is not star shaped");
- AllTris.erase(AllTris.begin());
- worst->forceRadius(-1);
- AllTris.insert(worst);
- delete v;
- }
- else
- gf->mesh_vertices.push_back(v);
- }
- else
- {
-// Msg(DEBUG,"Point %g %g is outside %g %g - %g %g - %g %g (%g %g)",
-// center[0],center[1],
-// Us [(worst)->tri()->getVertex(0)->getNum()],
-// Vs [(worst)->tri()->getVertex(0)->getNum()],
-// Us [(worst)->tri()->getVertex(1)->getNum()],
-// Vs [(worst)->tri()->getVertex(1)->getNum()],
-// Us [(worst)->tri()->getVertex(2)->getNum()],
-// Vs [(worst)->tri()->getVertex(2)->getNum()],
-// uv[0],uv[1]);
- AllTris.erase(AllTris.begin());
- worst->forceRadius(0);
- AllTris.insert(worst);
- // break;
- }
+ else {
+ // printf("%g %g %g\n",metric[0],metric[1],metric[2]);
+ // Msg(DEBUG,"Point %g %g is outside %g %g - %g %g - %g %g (%22.15E %22.15E)",
+ // center[0],center[1],
+ // Us [(worst)->tri()->getVertex(0)->getNum()],
+ // Vs [(worst)->tri()->getVertex(0)->getNum()],
+ // Us [(worst)->tri()->getVertex(1)->getNum()],
+ // Vs [(worst)->tri()->getVertex(1)->getNum()],
+ // Us [(worst)->tri()->getVertex(2)->getNum()],
+ // Vs [(worst)->tri()->getVertex(2)->getNum()],
+ // uv[0],uv[1]);
+ // throw;
+ AllTris.erase(AllTris.begin());
+ worst->forceRadius(0);
+ AllTris.insert(worst);
+ // break;
}
- }
-
+ }
+ }
+ for (int i=0;i<10;i++){
+ // bool splits = edgeSplitPass(1.4,gf,AllTris,SPCR_ALLWAYS,Us,Vs,vSizes,vSizesBGM);
+ // edgeSwapPass(gf,AllTris,SWCR_QUAL,Us,Vs,vSizes,vSizesBGM);
+ // bool collapses = edgeCollapsePass(0.67,gf,AllTris,Us,Vs,vSizes,vSizesBGM);
+ // edgeSwapPass(gf,AllTris,SWCR_QUAL,Us,Vs,vSizes,vSizesBGM);
+ // if (!collapses && !splits)break;
+ }
+
+ // for (int i=0;i<1;i++)if(!edgeCollapsePass(100,gf,AllTris,Us,Vs,vSizes,vSizesBGM))break;
+ // for (int i=0;i<100;i++)if(!edgeSwapPass(gf,AllTris,SWCR_QUAL,Us,Vs,vSizes,vSizesBGM))break;
+
+ // char name[245];
+ // sprintf(name,"paramFinal%d.pos",gf->tag());
+ // _printTris (name, AllTris, Us,Vs);
+ // _printTris ("yo.pos", AllTris, Us,Vs);
// optimize the mesh
-
+
// fill new gmsh structures with triangles
- while (1)
- {
- if (AllTris.begin() == AllTris.end() ) break;
- MTri3 *worst = *AllTris.begin();
- if (worst->isDeleted())
- delete worst->tri();
- else
- gf->triangles.push_back(worst->tri());
- delete worst;
- AllTris.erase(AllTris.begin());
- }
+ while (1) {
+ if (AllTris.begin() == AllTris.end() ) break;
+ MTri3 *worst = *AllTris.begin();
+ if (worst->isDeleted())
+ delete worst->tri();
+ else
+ gf->triangles.push_back(worst->tri());
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ }
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 5d307bc3f6a41e01144d2d16e60bea4cb6557bb8..9c10c3f49435670eac70913e434892cdd1ff49a3 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -30,19 +30,22 @@ class GFace;
class BDS_Mesh;
class BDS_Point;
-struct gmsh2dMetric
-{
- double a11,a21,a22;
- gmsh2dMetric (double lc); // uniform metric
- gmsh2dMetric (double _a11 = 1, double _a21= 0, double _a22=1)
- : a11(_a11),a21(_a21),a22(_a22)
- {}
- inline double eval ( const double &x, const double &y ) const
- {
- return x * a11 * x + 2 * x * a21 * y + y * a22 * y;
- }
-};
-
+void buildMetric ( GFace *gf , double *uv, double *metric);
+int inCircumCircleAniso ( GFace *gf, MTriangle *base, const double *uv , const double *metric,
+ const std::vector<double> & Us, const std::vector<double> & Vs);
+void circumCenterMetric ( MTriangle *base,
+ const double *metric,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs,
+ double *x, double &Radius2);
+bool circumCenterMetricInTriangle ( MTriangle *base,
+ const double *metric,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs);
+bool invMapUV ( MTriangle *t ,
+ double *p,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs , double *uv, double tol);
class MTri3
{
@@ -71,8 +74,6 @@ class MTri3
return inCircumCircle ( v->x(), v->y() );
}
- void Center_Circum_Aniso(double a, double b, double d, double &x, double &y, double &r) const ;
-
double getSurfaceXY () const { return base -> getSurfaceXY() ; };
double getSurfaceUV (GFace* gf) const { return base -> getSurfaceUV(gf) ; };
inline void setDeleted (bool d)
@@ -97,6 +98,7 @@ class MTri3
};
void connectTriangles ( std::list<MTri3*> & );
+void connectTriangles ( std::vector<MTri3*> & );
void insertVerticesInFace (GFace *gf, BDS_Mesh *);
class compareTri3Ptr
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 677b46dc83b8e2c9aa9d8c27ef6d560acf6a74f1..f460af3b7c1cef6bd5e11361c14987359cf37bf2 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1,19 +1,18 @@
#include "meshGFaceOptimize.h"
+#include "qualityMeasures.h"
#include "GFace.h"
#include "GEdge.h"
#include "GVertex.h"
#include "MVertex.h"
#include "MElement.h"
+#include "BackgroundMesh.h"
-typedef std::map<MVertex*,std::vector<MTriangle*> > v2t_cont ;
-typedef std::map<MEdge, std::pair<MTriangle*,MTriangle*> , Less_Edge> e2t_cont ;
-
-void buildVertexToTriangle ( GFace *gf , v2t_cont &adj )
+void buildVertexToTriangle ( std::vector<MTriangle*> &triangles, v2t_cont &adj )
{
adj.clear();
- for (unsigned int i=0;i<gf->triangles.size();i++)
+ for (unsigned int i=0;i<triangles.size();i++)
{
- MTriangle *t = gf->triangles[i];
+ MTriangle *t = triangles[i];
for (unsigned int j=0;j<3;j++)
{
MVertex *v = t->getVertex(j);
@@ -76,7 +75,7 @@ void parametricCoordinates ( MTriangle *t , GFace *gf, double u[3], double v[3])
void laplaceSmoothing (GFace *gf)
{
v2t_cont adj;
- buildVertexToTriangle ( gf , adj );
+ buildVertexToTriangle ( gf->triangles , adj );
for (int i=0;i<5;i++)
{
@@ -120,85 +119,553 @@ void laplaceSmoothing (GFace *gf)
extern void fourthPoint (double *p1, double *p2, double *p3, double *p4);
-bool edgeSwap(e2t_cont &adj, e2t_cont::iterator &it, GFace *gf)
-{
- MVertex *v1 = it->first.getVertex(0);
- MVertex *v2 = it->first.getVertex(1);
- MVertex *o1,*o2 ;
- MTriangle *t1 = it->second.first;
- MTriangle *t2 = it->second.second;
+double surfaceTriangleUV (MVertex *v1, MVertex *v2, MVertex *v3,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs){
+ const double v12[2] = {Us[v2->getNum()]-Us[v1->getNum()],Vs[v2->getNum()]-Vs[v1->getNum()]};
+ const double v13[2] = {Us[v3->getNum()]-Us[v1->getNum()],Vs[v3->getNum()]-Vs[v1->getNum()]};
+ return 0.5*fabs (v12[0]*v13[1]-v12[1]*v13[0]);
+}
+
+
+
+bool gmshEdgeSwap(MTri3 *t1,
+ GFace *gf,
+ int iLocalEdge,
+ std::vector<MTri3*> &newTris,
+ const gmshSwapCriterion &cr,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs,
+ const std::vector<double> & vSizes ,
+ const std::vector<double> & vSizesBGM ){
+
+ MTri3 *t2 = t1->getNeigh(iLocalEdge);
+ if (!t2) return false;
+
+ MVertex *v1 = t1->tri()->getVertex(iLocalEdge==0 ? 2 : iLocalEdge -1);
+ MVertex *v2 = t1->tri()->getVertex((iLocalEdge)%3);
+ MVertex *v3 = t1->tri()->getVertex((iLocalEdge+1)%3);
+ MVertex *v4 = 0 ;
for (int i=0;i<3;i++)
+ if (t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2)
+ v4 = t2->tri()->getVertex(i);
+
+ const double volumeRef = surfaceTriangleUV (v1,v2,v3,Us,Vs) + surfaceTriangleUV (v1,v2,v4,Us,Vs);
+
+ /// printf("%p %p %p %p\n",v2,v3,v4);
+ MTriangle *t1b = new MTriangle (v2,v3,v4);
+ /// printf("%p %p %p %p\n",v2,v3,v4);
+ MTriangle *t2b = new MTriangle (v4,v3,v1);
+
+ const double v1b = surfaceTriangleUV (v2,v3,v4,Us,Vs);
+ const double v2b = surfaceTriangleUV (v4,v3,v1,Us,Vs);
+ const double volume = v1b+v2b;
+ if (fabs(volume-volumeRef) > 1.e-10 * (volume+volumeRef) ||
+ v1b < 1.e-8 * (volume+volumeRef) ||
+ v2b < 1.e-8 * (volume+volumeRef)){
+ delete t1b;
+ delete t2b;
+ return false;
+ }
+
+ switch(cr){
+ case SWCR_QUAL:
{
- if (t1->getVertex(i) != v1 &&
- t1->getVertex(i) != v2)o1 = t1->getVertex(i);
- if (t2->getVertex(i) != v1 &&
- t2->getVertex(i) != v2)o1 = t2->getVertex(i);
+ const double triQualityRef = std::min(qmTriangle(t1->tri(),QMTRI_RHO),qmTriangle(t2->tri(),QMTRI_RHO));
+ const double triQuality = std::min(qmTriangle(t1b,QMTRI_RHO),qmTriangle(t2b,QMTRI_RHO));
+ if (triQuality < triQualityRef){
+ delete t1b;
+ delete t2b;
+ return false;
+ }
+ break;
}
-
+ case SWCR_DEL:
+ {
+ double edgeCenter[2] ={(Us[v1->getNum()]+Us[v2->getNum()])*.5,
+ (Vs[v1->getNum()]+Vs[v2->getNum()])*.5};
+ double uv4[2] ={Us[v4->getNum()],Vs[v4->getNum()]};
+ double metric[3];
+ buildMetric ( gf , edgeCenter , metric);
+ if (!inCircumCircleAniso (gf,t1->tri(),uv4,metric,Us,Vs)){
+ delete t1b;
+ delete t2b;
+ return false;
+ }
+ }
+ break;
+ case SWCR_CLOSE:
+ {
+ double avg1[3] = {(v1->x() + v2->x()) *.5,(v1->y() + v2->y()) *.5,(v1->z() + v2->z()) *.5};
+ double avg2[3] = {(v3->x() + v4->x()) *.5,(v3->y() + v4->y()) *.5,(v3->z() + v4->z()) *.5};
+
+ GPoint gp1 = gf->point(SPoint2((Us[v1->getNum()]+Us[v2->getNum()])*.5,(Vs[v1->getNum()]+Vs[v2->getNum()])*.5));
+ GPoint gp2 = gf->point(SPoint2((Us[v3->getNum()]+Us[v4->getNum()])*.5,(Vs[v3->getNum()]+Vs[v4->getNum()])*.5));
+ double d1 = (avg1[0]-gp1.x()) *(avg1[0]-gp1.x()) +(avg1[1]-gp1.y()) *(avg1[1]-gp1.y()) +(avg1[2]-gp1.z()) *(avg1[2]-gp1.z());
+ double d2 = (avg2[0]-gp2.x()) *(avg2[0]-gp2.x()) +(avg2[1]-gp2.y()) *(avg2[1]-gp2.y()) +(avg2[2]-gp2.z()) *(avg2[2]-gp2.z());
+ if (d1 < d2){
+ delete t1b;
+ delete t2b;
+ return false;
+ }
+ }
+ break;
+ default :
+ throw;
+ }
+
+ std::list<MTri3*> cavity;
+ for(int i=0;i<3;i++){
+ if (t1->getNeigh(i) && t1->getNeigh(i) != t2){
+ bool found = false;
+ for (std::list<MTri3*>::iterator it = cavity.begin();it!=cavity.end();it++){
+ if (*it == t1->getNeigh(i))found = true;
+ }
+ if (!found)cavity.push_back(t1->getNeigh(i));
+ }
+ }
+ for(int i=0;i<3;i++){
+ if (t2->getNeigh(i) && t2->getNeigh(i) != t1){
+ bool found = false;
+ for (std::list<MTri3*>::iterator it = cavity.begin();it!=cavity.end();it++){
+ if (*it == t2->getNeigh(i))found = true;
+ }
+ if (!found)cavity.push_back(t2->getNeigh(i));
+ }
+ }
+ double lc1 = 0.3333333333*(vSizes [t1b->getVertex(0)->getNum()]+
+ vSizes [t1b->getVertex(1)->getNum()]+
+ vSizes [t1b->getVertex(2)->getNum()]);
+ double lcBGM1 = 0.3333333333*(vSizesBGM [t1b->getVertex(0)->getNum()]+
+ vSizesBGM [t1b->getVertex(1)->getNum()]+
+ vSizesBGM [t1b->getVertex(2)->getNum()]);
+ double lc2 = 0.3333333333*(vSizes [t2b->getVertex(0)->getNum()]+
+ vSizes [t2b->getVertex(1)->getNum()]+
+ vSizes [t2b->getVertex(2)->getNum()]);
+ double lcBGM2 = 0.3333333333*(vSizesBGM [t2b->getVertex(0)->getNum()]+
+ vSizesBGM [t2b->getVertex(1)->getNum()]+
+ vSizesBGM [t2b->getVertex(2)->getNum()]);
+ MTri3 *t1b3 = new MTri3 ( t1b , Extend1dMeshIn2dSurfaces() ? std::min(lc1,lcBGM1) : lcBGM1 );
+ MTri3 *t2b3 = new MTri3 ( t2b , Extend1dMeshIn2dSurfaces() ? std::min(lc2,lcBGM2) : lcBGM2 );
+
+ cavity.push_back(t1b3);
+ cavity.push_back(t2b3);
+ t1->setDeleted(true);
+ t2->setDeleted(true);
+ connectTriangles ( cavity );
+ newTris.push_back(t1b3);
+ newTris.push_back(t2b3);
+ return true;
}
+inline double computeEdgeAdimLength(MVertex *v1, MVertex *v2, GFace *f,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs,
+ const std::vector<double> & vSizes ,
+ const std::vector<double> & vSizesBGM ){
+ const double edgeCenter[2] ={(Us[v1->getNum()]+Us[v2->getNum()])*.5,
+ (Vs[v1->getNum()]+Vs[v2->getNum()])*.5};
+ GPoint GP = f->point (edgeCenter[0],edgeCenter[1]);
+
+ const double dx1 = v1->x()-GP.x();
+ const double dy1 = v1->y()-GP.y();
+ const double dz1 = v1->z()-GP.z();
+ const double l1 = sqrt(dx1*dx1+dy1*dy1+dz1*dz1);
+ const double dx2 = v2->x()-GP.x();
+ const double dy2 = v2->y()-GP.y();
+ const double dz2 = v2->z()-GP.z();
+ const double l2 = sqrt(dx2*dx2+dy2*dy2+dz2*dz2);
+ if (Extend1dMeshIn2dSurfaces())return 2*(l1+l2)/(std::min(vSizes[v1->getNum()],vSizesBGM[v1->getNum()]) +
+ std::min(vSizes[v2->getNum()],vSizesBGM[v2->getNum()]));
+ return 2*(l1+l2)/(vSizesBGM[v1->getNum()] + vSizesBGM[v2->getNum()]);
+}
+bool gmshEdgeSplit(const double lMax,
+ MTri3 *t1,
+ GFace *gf,
+ int iLocalEdge,
+ std::vector<MTri3*> &newTris,
+ const gmshSplitCriterion &cr,
+ std::vector<double> & Us,
+ std::vector<double> & Vs,
+ std::vector<double> & vSizes ,
+ std::vector<double> & vSizesBGM ){
+
+ MTri3 *t2 = t1->getNeigh(iLocalEdge);
+ if (!t2) return false;
-bool edgeSwapTestDelaunay(e2t_cont::iterator &it,GFace *gf)
-{
- if(!it->second.second) return false;
+ MVertex *v1 = t1->tri()->getVertex(iLocalEdge==0 ? 2 : iLocalEdge -1);
+ MVertex *v2 = t1->tri()->getVertex((iLocalEdge)%3);
+ double edgeCenter[2] ={(Us[v1->getNum()]+Us[v2->getNum()])*.5,
+ (Vs[v1->getNum()]+Vs[v2->getNum()])*.5};
+
+ double al = computeEdgeAdimLength(v1, v2, gf,Us,Vs,vSizes,vSizesBGM);
+ if (al <= lMax )return false;
+ // printf("al,lMax %12.5E %12.5E \n",al,lMax);
- MVertex *v1 = it->first.getVertex(0);
- MVertex *v2 = it->first.getVertex(1);
- MVertex *o1,*o2 ;
+ MVertex *v3 = t1->tri()->getVertex((iLocalEdge+1)%3);
+ MVertex *v4 = 0 ;
for (int i=0;i<3;i++)
+ if (t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2)
+ v4 = t2->tri()->getVertex(i);
+
+ GPoint p = gf->point (edgeCenter[0],edgeCenter[1]);
+ MVertex *vnew = new MFaceVertex (p.x(),p.y(),p.z(),gf,edgeCenter[0],edgeCenter[1]);
+
+ MTriangle *t1b = new MTriangle (v1,vnew,v3);
+ MTriangle *t2b = new MTriangle (vnew,v2,v3);
+ MTriangle *t3b = new MTriangle (v1,v4,vnew);
+ MTriangle *t4b = new MTriangle (v4,v2,vnew);
+
+ switch(cr){
+ case SPCR_QUAL:
{
- if (it->second.first->getVertex(i) != v1 &&
- it->second.first->getVertex(i) != v2)o1 = it->second.first->getVertex(i);
- if (it->second.second->getVertex(i) != v1 &&
- it->second.second->getVertex(i) != v2)o2 = it->second.second->getVertex(i);
- }
- double p1[2];
- double p2[2];
- double op1[2];
- double op2[2];
- parametricCoordinates ( v1,gf,p1[0],p1[1] );
- parametricCoordinates ( v2,gf,p2[0],p2[1] );
- parametricCoordinates ( o1,gf,op1[0],op1[1] );
- parametricCoordinates ( o2,gf,op2[0],op2[1] );
+ const double triQualityRef = std::min(qmTriangle(t1->tri(),QMTRI_RHO),qmTriangle(t2->tri(),QMTRI_RHO));
+ const double triQuality =
+ std::min(std::min(std::min(qmTriangle(t1b,QMTRI_RHO),qmTriangle(t2b,QMTRI_RHO)),
+ qmTriangle(t3b,QMTRI_RHO)),qmTriangle(t4b,QMTRI_RHO));
+ if (triQuality < triQualityRef){
+ delete t1b;
+ delete t2b;
+ delete t3b;
+ delete t4b;
+ delete vnew;
+ return false;
+ }
+ break;
+ }
+ case SPCR_ALLWAYS:
+ break;
+ default :
+ throw;
+ }
+
+ gf->mesh_vertices.push_back(vnew);
+ vnew->setNum(Us.size());
+ double lcN = 0.5 * (vSizes [v1->getNum()] + vSizes [v2->getNum()] );
+ double lcBGM = BGM_MeshSize(gf,edgeCenter[0],edgeCenter[1],p.x(),p.y(),p.z());
+
+ vSizesBGM.push_back( lcBGM );
+ vSizes.push_back ( lcN);
+ Us.push_back( edgeCenter[0] );
+ Vs.push_back( edgeCenter[1] );
+
+ std::list<MTri3*> cavity;
+ for(int i=0;i<3;i++){
+ if (t1->getNeigh(i) && t1->getNeigh(i) != t2){
+ bool found = false;
+ for (std::list<MTri3*>::iterator it = cavity.begin();it!=cavity.end();it++){
+ if (*it == t1->getNeigh(i))found = true;
+ }
+ if (!found)cavity.push_back(t1->getNeigh(i));
+ }
+ }
+ for(int i=0;i<3;i++){
+ if (t2->getNeigh(i) && t2->getNeigh(i) != t1){
+ bool found = false;
+ for (std::list<MTri3*>::iterator it = cavity.begin();it!=cavity.end();it++){
+ if (*it == t2->getNeigh(i))found = true;
+ }
+ if (!found)cavity.push_back(t2->getNeigh(i));
+ }
+ }
+ double lc1 = 0.3333333333*(vSizes [t1b->getVertex(0)->getNum()]+
+ vSizes [t1b->getVertex(1)->getNum()]+
+ vSizes [t1b->getVertex(2)->getNum()]);
+ double lcBGM1 = 0.3333333333*(vSizesBGM [t1b->getVertex(0)->getNum()]+
+ vSizesBGM [t1b->getVertex(1)->getNum()]+
+ vSizesBGM [t1b->getVertex(2)->getNum()]);
+ double lc2 = 0.3333333333*(vSizes [t2b->getVertex(0)->getNum()]+
+ vSizes [t2b->getVertex(1)->getNum()]+
+ vSizes [t2b->getVertex(2)->getNum()]);
+ double lcBGM2 = 0.3333333333*(vSizesBGM [t2b->getVertex(0)->getNum()]+
+ vSizesBGM [t2b->getVertex(1)->getNum()]+
+ vSizesBGM [t2b->getVertex(2)->getNum()]);
+ double lc3 = 0.3333333333*(vSizes [t3b->getVertex(0)->getNum()]+
+ vSizes [t3b->getVertex(1)->getNum()]+
+ vSizes [t3b->getVertex(2)->getNum()]);
+ double lcBGM3 = 0.3333333333*(vSizesBGM [t3b->getVertex(0)->getNum()]+
+ vSizesBGM [t3b->getVertex(1)->getNum()]+
+ vSizesBGM [t3b->getVertex(2)->getNum()]);
+ double lc4 = 0.3333333333*(vSizes [t4b->getVertex(0)->getNum()]+
+ vSizes [t4b->getVertex(1)->getNum()]+
+ vSizes [t4b->getVertex(2)->getNum()]);
+ double lcBGM4 = 0.3333333333*(vSizesBGM [t4b->getVertex(0)->getNum()]+
+ vSizesBGM [t4b->getVertex(1)->getNum()]+
+ vSizesBGM [t4b->getVertex(2)->getNum()]);
+ MTri3 *t1b3 = new MTri3 ( t1b , Extend1dMeshIn2dSurfaces() ? std::min(lc1,lcBGM1) : lcBGM1 );
+ MTri3 *t2b3 = new MTri3 ( t2b , Extend1dMeshIn2dSurfaces() ? std::min(lc2,lcBGM2) : lcBGM2 );
+ MTri3 *t3b3 = new MTri3 ( t3b , Extend1dMeshIn2dSurfaces() ? std::min(lc3,lcBGM3) : lcBGM3 );
+ MTri3 *t4b3 = new MTri3 ( t4b , Extend1dMeshIn2dSurfaces() ? std::min(lc4,lcBGM4) : lcBGM4 );
+
+ cavity.push_back(t1b3);
+ cavity.push_back(t2b3);
+ cavity.push_back(t3b3);
+ cavity.push_back(t4b3);
+ t1->setDeleted(true);
+ t2->setDeleted(true);
+ connectTriangles ( cavity );
+ newTris.push_back(t1b3);
+ newTris.push_back(t2b3);
+ newTris.push_back(t3b3);
+ newTris.push_back(t4b3);
+
+ return true;
+}
+
+void computeNeighboringTrisOfACavity (const std::vector<MTri3*> &cavity,std::vector<MTri3*> &outside)
+{
+ outside.clear();
+ for (int i=0;i<cavity.size();i++){
+ for (int j=0;j<3;j++){
+ MTri3 * neigh = cavity[i]->getNeigh(j);
+ if (neigh)
+ {
+ bool found = false;
+ for (int k=0;k<outside.size();k++){
+ if(outside[k]==neigh){
+ found=true;
+ break;
+ }
+ }
+ if (!found){
+ for (int k=0;k<cavity.size() ;k++){
+ if(cavity[k] ==neigh){
+ found=true;
+ }
+ }
+ }
+ if(!found)outside.push_back(neigh);
+ }
+ }
+ }
+}
+bool gmshBuildVertexCavity ( MTri3 *t,
+ int iLocalVertex,
+ MVertex **v1,
+ std::vector<MTri3*> &cavity,
+ std::vector<MTri3*> &outside,
+ std::vector<MVertex*> &ring ){
+ cavity.clear();
+ ring.clear();
+
+ *v1 = t->tri()->getVertex(iLocalVertex);
+
+// printf("VERTEX %d\n",
+// t->tri()->getVertex(iLocalVertex)->getNum());
+
+ MVertex *lastinring = t->tri()->getVertex((iLocalVertex+1)%3);
+ ring.push_back(lastinring);
+ cavity.push_back(t);
+
+// printf("START triangle %d %d %d, vertex %d\n",
+// t->tri()->getVertex(0)->getNum(),
+// t->tri()->getVertex(1)->getNum(),
+// t->tri()->getVertex(2)->getNum(),
+// lastinring->getNum());
+
+
+ while (1){
+ int iEdge = -1;
+ // printf("look for %d %d\n",(*v1)->getNum(),lastinring->getNum());
+ for (int i=0;i<3;i++){
+ MVertex *v2 = t->tri()->getVertex((i+2)%3);
+ MVertex *v3 = t->tri()->getVertex(i);
+ // printf("--> %d %d\n",v2->getNum(),v3->getNum());
+ if ( (v2 == *v1 && v3 == lastinring ) ||
+ (v2 == lastinring && v3 == *v1 )){
+ iEdge = i;
+ t = t->getNeigh(i);
+ if (t==cavity[0]) {
+ computeNeighboringTrisOfACavity (cavity,outside);
+ return true;
+ }
+ if (!t)return false;
+ if (t->isDeleted()){printf("weird\n");throw;}
+ cavity.push_back(t);
+ for (int j=0;j<3;j++){
+ if (t->tri()->getVertex(j) !=lastinring && t->tri()->getVertex(j) != *v1){
+ lastinring = t->tri()->getVertex(j);
+ ring.push_back(lastinring);
+ j = 100;
+ }
+ }
+// printf("CONTINUE (%d) triangle %p = %d %d %d, vertex %d size %d\n",i,
+// t,
+// t->tri()->getVertex(0)->getNum(),
+// t->tri()->getVertex(1)->getNum(),
+// t->tri()->getVertex(2)->getNum(),
+// lastinring->getNum(),ring.size());
+ break;
+ }
+ }
+ if (iEdge == -1) {printf("not found\n"); throw;}
+ }
+}
+
+
+
+bool gmshVertexCollapse(const double lMin,
+ MTri3 *t1,
+ GFace *gf,
+ int iLocalVertex,
+ std::vector<MTri3*> &newTris,
+ std::vector<double> & Us,
+ std::vector<double> & Vs,
+ std::vector<double> & vSizes ,
+ std::vector<double> & vSizesBGM ){
+ MVertex *v;
+ std::vector<MTri3*> cavity;
+ std::vector<MTri3*> outside;
+ std::vector<MVertex*> ring ;
+ // printf("%p \n",t1);
+ if (!gmshBuildVertexCavity (t1,iLocalVertex, &v,cavity,outside,ring) )return false;
- double ori_t1 = gmsh::orient2d(op1, p1, op2);
- double ori_t2 = gmsh::orient2d(op1,op2, p2);
+ double l_min = lMin;
+ int iMin = -1;
+ for (int i=0;i<ring.size();i++){
+ double l = computeEdgeAdimLength(v, ring[i], gf,Us,Vs,vSizes,vSizesBGM);
+ if (l < l_min){
+ iMin = i;
+ }
+ }
+ if (iMin == -1) return false;
+
+ double surfBefore = 0.0;
+ for (int i=0;i<ring.size();i++){
+ MVertex *v1 = ring[i];
+ MVertex *v2 = ring[(i + 1)%ring.size()];
+ surfBefore +=surfaceTriangleUV (v1,v2,v,Us,Vs);
+ }
+ double surfAfter = 0.0;
+ for (int i=0;i<ring.size()-2;i++){
+ MVertex *v1 = ring[(iMin + 1 + i)%ring.size()];
+ MVertex *v2 = ring[(iMin + 1 + i + 1)%ring.size()];
+ double sAfter = surfaceTriangleUV (v1,v2,ring[iMin],Us,Vs);
+ double sBefore = surfaceTriangleUV (v1,v2,v,Us,Vs);
+ if (sAfter < 0.1 * sBefore)return false;
+ surfAfter += sAfter;
+ }
- // the quad is concave in parametric coord, return false
- if ( ori_t1 * ori_t2 <= 0) return false;
+ // printf("%12.5E %12.5E %d\n",surfBefore,surfAfter,iMin);
+ if (fabs(surfBefore - surfAfter) > 1.e-10*(surfBefore + surfAfter))return false;
- double p1x[3] = {v1->x(),v1->y(),v1->z()};
- double p2x[3] = {v2->x(),v2->y(),v2->z()};
- double op1x[3] = {o1->x(),o1->y(),o1->z()};
- double op2x[3] = {o2->x(),o2->y(),o2->z()};
- double fourth[3];
- fourthPoint(p1x,p2x,op1x,fourth);
- double result = gmsh::insphere(p1x, p2x, op1x, fourth, op2x) * gmsh::orient3d(p1x, p2x, op1x, fourth);
- return result > 0.;
+ for (int i=0;i<ring.size()-2;i++){
+ MVertex *v1 = ring[(iMin + 1 + i)%ring.size()];
+ MVertex *v2 = ring[(iMin + 1 + i + 1)%ring.size()];
+ MTriangle *t = new MTriangle(v1,v2,ring[iMin]);
+ double lc = 0.3333333333*(vSizes [t->getVertex(0)->getNum()]+
+ vSizes [t->getVertex(1)->getNum()]+
+ vSizes [t->getVertex(2)->getNum()]);
+ double lcBGM = 0.3333333333*(vSizesBGM [t->getVertex(0)->getNum()]+
+ vSizesBGM [t->getVertex(1)->getNum()]+
+ vSizesBGM [t->getVertex(2)->getNum()]);
+ MTri3 *t3 = new MTri3 ( t , Extend1dMeshIn2dSurfaces() ? std::min(lc,lcBGM) : lcBGM);
+ // printf("Creation %p = %d %d %d\n",t3,v1->getNum(),v2->getNum(),ring[iMin]->getNum());
+ outside.push_back(t3);
+ newTris.push_back(t3);
+ }
+ for (int i=0;i<cavity.size();i++)
+ cavity[i]->setDeleted(true);
+ connectTriangles ( outside );
+ return true;
}
-void edgeSwappingLawson (GFace *gf)
+int edgeSwapPass (GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
+ const gmshSwapCriterion &cr,
+ const std::vector<double> & Us ,
+ const std::vector<double> & Vs,
+ const std::vector<double> & vSizes ,
+ const std::vector<double> & vSizesBGM)
{
+ typedef std::set<MTri3*,compareTri3Ptr> CONTAINER ;
+ std::vector<MTri3*> newTris;
+
+ int nbSwap = 0;
+
+ for (CONTAINER::iterator it = allTris.begin();it!=allTris.end();++it){
+ if (!(*it)->isDeleted()){
+ for (int i=0;i<3;i++){
+ if (gmshEdgeSwap(*it,gf,i,newTris,cr,Us,Vs,vSizes,vSizesBGM)) {nbSwap++;break;}
+ }
+ }
+ else{
+ CONTAINER::iterator itb = it;
+ ++it;
+ delete *itb;
+ allTris.erase(itb);
+ }
+ }
+ // printf("B %d %d tris ",allTris.size(),newTris.size());
+ allTris.insert(newTris.begin(),newTris.end());
+ // printf("A %d %d tris\n",allTris.size(),newTris.size());
+ return nbSwap;
+}
- e2t_cont adj;
- buildEdgeToTriangle ( gf , adj );
+int edgeSplitPass (double maxLC,
+ GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
+ const gmshSplitCriterion &cr,
+ std::vector<double> & Us ,
+ std::vector<double> & Vs,
+ std::vector<double> & vSizes ,
+ std::vector<double> & vSizesBGM)
+{
+ typedef std::set<MTri3*,compareTri3Ptr> CONTAINER ;
+ std::vector<MTri3*> newTris;
- e2t_cont :: iterator it = adj.begin();
+ int nbSplit = 0;
- while (it != adj.end())
- {
- bool evalSwapLawson = edgeSwapTestDelaunay(it,gf);
- if ( evalSwapLawson )
- {
- edgeSwap (adj, it ,gf );
- e2t_cont :: iterator itb = it;
- ++it;
- adj.erase(itb);
- }
- else
- ++it;
- }
+ for (CONTAINER::iterator it = allTris.begin();it!=allTris.end();++it){
+ if (!(*it)->isDeleted()){
+ for (int i=0;i<3;i++){
+ if (gmshEdgeSplit (maxLC,*it,gf,i,newTris,cr,Us,Vs,vSizes,vSizesBGM)) {nbSplit++;break;}
+ }
+ }
+ else{
+ CONTAINER::iterator itb = it;
+ delete *it;
+ ++it;
+ allTris.erase(itb);
+ if (it == allTris.end())break;
+ }
+ }
+ printf("B %d %d tris ",allTris.size(),newTris.size());
+ allTris.insert(newTris.begin(),newTris.end());
+ printf("A %d %d tris\n",allTris.size(),newTris.size());
+ return nbSplit;
+}
+
+int edgeCollapsePass (double minLC,
+ GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
+ std::vector<double> & Us ,
+ std::vector<double> & Vs,
+ std::vector<double> & vSizes ,
+ std::vector<double> & vSizesBGM)
+{
+ typedef std::set<MTri3*,compareTri3Ptr> CONTAINER ;
+ std::vector<MTri3*> newTris;
+
+ int nbCollapse = 0;
+
+ for (CONTAINER::reverse_iterator it = allTris.rbegin();it!=allTris.rend();++it){
+ if (!(*it)->isDeleted()){
+ for (int i=0;i<3;i++){
+ if (gmshVertexCollapse (minLC,*it,gf,i,newTris,Us,Vs,vSizes,vSizesBGM)) {nbCollapse++; break;}
+ }
+ }
+// else{
+// CONTAINER::reverse_iterator itb = it;
+// delete *it;
+// ++it;
+// allTris.rerase(itb);
+// if (it == allTris.rend())break;
+// }
+// if (nbCollapse == 114)break;
+ }
+ printf("B %d %d tris ",allTris.size(),newTris.size());
+ allTris.insert(newTris.begin(),newTris.end());
+ printf("A %d %d tris\n",allTris.size(),newTris.size());
+ return nbCollapse;
}
diff --git a/Mesh/meshGFaceOptimize.h b/Mesh/meshGFaceOptimize.h
index b68331f14c015fea07265b220d4ebc4c271dac6f..e3f0d8e404178268018f45735afd0383f982380b 100644
--- a/Mesh/meshGFaceOptimize.h
+++ b/Mesh/meshGFaceOptimize.h
@@ -21,6 +21,7 @@
// Please report all bugs and problems to <gmsh@geuz.org>.
#include "MElement.h"
#include "MEdge.h"
+#include "meshGFaceDelaunayInsertion.h"
#include <map>
#include <vector>
class GFace;
@@ -31,4 +32,49 @@ void buildVertexToTriangle ( std::vector<MTriangle*> & , v2t_cont &adj );
void buildEdgeToTriangle ( std::vector<MTriangle*> & , e2t_cont &adj );
void laplaceSmoothing (GFace *gf);
void edgeSwappingLawson (GFace *gf);
+enum gmshSwapCriterion {SWCR_DEL, SWCR_QUAL,SWCR_NORM,SWCR_CLOSE};
+enum gmshSplitCriterion {SPCR_CLOSE, SPCR_QUAL,SPCR_ALLWAYS};
+int edgeSwapPass (GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
+ const gmshSwapCriterion &cr,
+ const std::vector<double> & Us ,
+ const std::vector<double> & Vs,
+ const std::vector<double> & vSizes ,
+ const std::vector<double> & vSizesBGM);
+
+
+bool gmshEdgeSwap(MTri3 *t1,
+ GFace *gf,
+ int iLocalEdge,
+ std::vector<MTri3*> &newTris,
+ const gmshSwapCriterion &cr,
+ const std::vector<double> & Us,
+ const std::vector<double> & Vs,
+ const std::vector<double> & vSizes,
+ const std::vector<double> & vSizesBGM);
+
+bool gmshEdgeSplit(const double lMax,
+ MTri3 *t1,
+ GFace *gf,
+ int iLocalEdge,
+ std::vector<MTri3*> &newTris,
+ const gmshSplitCriterion &cr,
+ std::vector<double> & Us,
+ std::vector<double> & Vs,
+ std::vector<double> & vSizes,
+ std::vector<double> & vSizesBGM);
+
+int edgeSplitPass (double maxLC,
+ GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
+ const gmshSplitCriterion &cr,
+ std::vector<double> & Us ,
+ std::vector<double> & Vs,
+ std::vector<double> & vSizes ,
+ std::vector<double> & vSizesBGM);
+
+int edgeCollapsePass (double minLC,
+ GFace *gf, std::set<MTri3*,compareTri3Ptr> &allTris,
+ std::vector<double> & Us ,
+ std::vector<double> & Vs,
+ std::vector<double> & vSizes ,
+ std::vector<double> & vSizesBGM);
#endif
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 0635f2abdd3b736672ff8da52ee3708dd610d49a..6cd56c4c20ec6ddbc0f2d9dc518ec47688148359 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.37 2007-11-28 14:18:10 remacle Exp $
+// $Id: meshGRegion.cpp,v 1.38 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -609,3 +609,73 @@ void optimizeMeshGRegionGmsh::operator() (GRegion *gr)
gmshOptimizeMesh (gr, QMTET_2);
}
+
+bool buildFaceSearchStructure ( GModel *model , fs_cont&search )
+{
+ search.clear();
+
+ GModel::fiter fit = model->firstFace() ;
+ while (fit != model->lastFace()){
+ for (int i=0;i<(*fit)->triangles.size();i++)
+ {
+ MVertex *p1=(*fit)->triangles[i]->getVertex(0);
+ MVertex *p2=(*fit)->triangles[i]->getVertex(1);
+ MVertex *p3=(*fit)->triangles[i]->getVertex(2);
+ MVertex *p = std::min(p1,std::min(p2,p3));
+ search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p, std::pair<MTriangle*,GFace*>((*fit)->triangles[i],*fit)));
+ }
+ ++fit;
+ }
+ return true;
+}
+
+bool buildEdgeSearchStructure ( GModel *model , es_cont&search )
+{
+ search.clear();
+
+ GModel::eiter eit = model->firstEdge() ;
+ while (eit != model->lastEdge()){
+ for (int i=0;i<(*eit)->lines.size();i++)
+ {
+ MVertex *p1=(*eit)->lines[i]->getVertex(0);
+ MVertex *p2=(*eit)->lines[i]->getVertex(1);
+ MVertex *p = std::min(p1,p2);
+ search.insert ( std::pair<MVertex*,std::pair<MLine*,GEdge*> > ( p, std::pair<MLine*,GEdge*>((*eit)->lines[i],*eit)));
+ }
+ ++eit;
+ }
+ return true;
+}
+
+GFace* findInFaceSearchStructure ( MVertex *p1,MVertex *p2,MVertex *p3, const fs_cont&search )
+{
+ MVertex *p = std::min(p1,std::min(p2,p3));
+
+ for (fs_cont::const_iterator it = search.lower_bound(p);
+ it != search.upper_bound(p);
+ ++it)
+ {
+ MTriangle *t = it->second.first;
+ GFace *gf = it->second.second;
+ if ((t->getVertex(0) == p1 ||t->getVertex(0) == p2 ||t->getVertex(0) == p3)&&
+ (t->getVertex(1) == p1 ||t->getVertex(1) == p2 ||t->getVertex(1) == p3)&&
+ (t->getVertex(2) == p1 ||t->getVertex(2) == p2 ||t->getVertex(2) == p3))return gf;
+ }
+ return 0;
+}
+
+GEdge* findInEdgeSearchStructure ( MVertex *p1, MVertex *p2, const es_cont&search )
+{
+ MVertex *p = std::min(p1,p2);
+
+ for (es_cont::const_iterator it = search.lower_bound(p);
+ it != search.upper_bound(p);
+ ++it)
+ {
+ MLine *l = it->second.first;
+ GEdge *ge = it->second.second;
+ if ((l->getVertex(0) == p1 ||l->getVertex(0) == p2)&&
+ (l->getVertex(1) == p1 ||l->getVertex(1) == p2))return ge;
+ }
+ return 0;
+}
diff --git a/Mesh/meshGRegion.h b/Mesh/meshGRegion.h
index 6f3fb49c7bfe205f396ea61558d6f492431d1237..f627658e648041e4228a79807270c47318fe10cc 100644
--- a/Mesh/meshGRegion.h
+++ b/Mesh/meshGRegion.h
@@ -21,10 +21,15 @@
// Please report all bugs and problems to <gmsh@geuz.org>.
#include <vector>
+#include <map>
class GModel;
class GRegion;
-
+class GFace;
+class GEdge;
+class MVertex;
+class MLine;
+class MTriangle;
// Create the mesh of the region
class meshGRegion {
public :
@@ -56,16 +61,24 @@ class deMeshGRegion {
void operator () (GRegion *);
};
-// adapt the mesh of a region
-class adaptMeshGRegion {
- public :
- void operator () (GRegion *);
-};
-
void MeshDelaunayVolume(std::vector<GRegion*> &delaunay);
int MeshTransfiniteVolume(GRegion *gr);
int SubdivideExtrudedMesh(GModel *m);
void carveHole(GRegion *gr, int num, double distance, std::vector<int> &surfaces);
+typedef std::multimap<MVertex*,std::pair<MTriangle*,GFace*> > fs_cont ;
+typedef std::multimap<MVertex*,std::pair<MLine*,GEdge*> > es_cont ;
+GFace* findInFaceSearchStructure ( MVertex *p1,MVertex *p2,MVertex *p3, const fs_cont&search );
+GEdge* findInEdgeSearchStructure ( MVertex *p1,MVertex *p2, const es_cont&search );
+bool buildFaceSearchStructure ( GModel *model , fs_cont&search );
+bool buildEdgeSearchStructure ( GModel *model , es_cont&search );
+
+// adapt the mesh of a region
+class adaptMeshGRegion {
+ es_cont *_es; fs_cont *_fs;
+ public :
+ void operator () (GRegion *);
+};
+
#endif
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 42dffd8cad07a94903a749bd4fe40a6ede67a400..d78c2e2104c7b07aff09072c8d418dd0fd8a7bcc 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.25 2007-12-03 15:17:40 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.26 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -21,11 +21,11 @@
#include "OS.h"
#include "BackgroundMesh.h"
+#include "meshGRegion.h"
#include "meshGRegionLocalMeshMod.h"
#include "meshGRegionDelaunayInsertion.h"
#include "GModel.h"
#include "GRegion.h"
-#include "meshGRegion.h"
#include "Numeric.h"
#include "Message.h"
#include <set>
@@ -327,42 +327,7 @@ static void setLcs ( MTetrahedron *t, std::map<MVertex*,double> &vSizes)
-GFace* findInFaceSearchStructure ( MVertex *p1,MVertex *p2,MVertex *p3, const fs_cont&search )
-{
- MVertex *p = std::min(p1,std::min(p2,p3));
-
- for (fs_cont::const_iterator it = search.lower_bound(p);
- it != search.upper_bound(p);
- ++it)
- {
- MTriangle *t = it->second.first;
- GFace *gf = it->second.second;
- if ((t->getVertex(0) == p1 ||t->getVertex(0) == p2 ||t->getVertex(0) == p3)&&
- (t->getVertex(1) == p1 ||t->getVertex(1) == p2 ||t->getVertex(1) == p3)&&
- (t->getVertex(2) == p1 ||t->getVertex(2) == p2 ||t->getVertex(2) == p3))return gf;
- }
- return 0;
-}
-
-bool buildFaceSearchStructure ( GModel *model , fs_cont&search )
-{
- search.clear();
-
- GModel::fiter fit = model->firstFace() ;
- while (fit != model->lastFace()){
- for (int i=0;i<(*fit)->triangles.size();i++)
- {
- MVertex *p1=(*fit)->triangles[i]->getVertex(0);
- MVertex *p2=(*fit)->triangles[i]->getVertex(1);
- MVertex *p3=(*fit)->triangles[i]->getVertex(2);
- MVertex *p = std::min(p1,std::min(p2,p3));
- search.insert ( std::pair<MVertex*,std::pair<MTriangle*,GFace*> > ( p, std::pair<MTriangle*,GFace*>((*fit)->triangles[i],*fit)));
- }
- ++fit;
- }
- return true;
-}
GRegion *getRegionFromBoundingFaces (GModel *model,
std::set<GFace *> &faces_bound)
@@ -471,7 +436,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
}
}
- Msg(INFO,"Opti : START with %12.5E QBAD %12.5E QAVG %12.5E",totalVolumeb,worst,avg/count);
+ Msg(INFO,"Adaptation : START with %12.5E QBAD %12.5E QAVG %12.5E",totalVolumeb,worst,avg/count);
for (int i=0;i<nbRanges;i++){
double low = (double)i/(nbRanges);
double high = (double)(i+1)/(nbRanges);
@@ -614,7 +579,9 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
typedef std::list<MTet4 *> CONTAINER ;
CONTAINER allTets;
for (int i=0;i<gr->tetrahedra.size();i++){
- allTets.push_back(new MTet4(gr->tetrahedra[i],qm));
+ MTet4 * t = new MTet4(gr->tetrahedra[i],qm);
+ t->setOnWhat(gr);
+ allTets.push_back(t);
}
gr->tetrahedra.clear();
@@ -655,7 +622,7 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
}
double qMin = 0.5;
- double sliverLimit = 0.2;
+ double sliverLimit = 0.1;
int nbESwap=0, nbFSwap=0, nbReloc=0;
@@ -695,7 +662,21 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
}
}
- if (!newTets.size())break;
+ if (!newTets.size()){
+ int nbSlivers = 0;
+ int nbSliversWeCanDoSomething = 0;
+ for (int i=0;i<illegals.size();i++)
+ if (!(illegals[i]->isDeleted())){
+ if (gmshSliverRemoval(newTets,illegals[i],qm))
+ nbSliversWeCanDoSomething++;
+ nbSlivers ++;
+ }
+ Msg(INFO,"Opti : %d Sliver Removals",nbSliversWeCanDoSomething);
+ }
+
+ if (!newTets.size()){
+ break;
+ }
// add all the new tets in the container
for (int i=0;i<newTets.size();i++){
@@ -739,15 +720,12 @@ void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm)
Msg(INFO,"Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",nbESwap,nbFSwap,nbReloc,totalVolumeb,worst,avg/count,t2-t1);
}
- int nbSlivers = 0;
- for (int i=0;i<illegals.size();i++)
- if (!(illegals[i]->isDeleted()))nbSlivers ++;
-
- if (nbSlivers){
- Msg(INFO,"Opti : %d illegal tets are still in the mesh, trying to remove them",nbSlivers);
+
+ if (illegals.size()){
+ Msg(INFO,"Opti : %d illegal tets are still in the mesh",illegals.size());
}
else{
- Msg(INFO,"Opti : no illegal tets in the mesh ;-)",nbSlivers);
+ Msg(INFO,"Opti : no illegal tets in the mesh ;-)");
}
for (int i=0;i<nbRanges;i++){
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 1c6469b5357e188f9abe373e21115bdaa225fdcf..30c002b89507d8a4a36f2d7f7fd113984a8f0643 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -235,8 +235,5 @@ private:
};
void gmshOptimizeMesh (GRegion *gr, const gmshQualityMeasure4Tet &qm);
-typedef std::multimap<MVertex*,std::pair<MTriangle*,GFace*> > fs_cont ;
-GFace* findInFaceSearchStructure ( MVertex *p1,MVertex *p2,MVertex *p3, const fs_cont&search );
-bool buildFaceSearchStructure ( GModel *model , fs_cont&search );
#endif
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index c2066916a715dcfac70098e979be01c92457f982..a3d4a919432c8c7f46c4f11ddf839dd087f25d25 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -1,5 +1,6 @@
#include "meshGRegionLocalMeshMod.h"
#include "GEntity.h"
+#include "GRegion.h"
#include "Message.h"
static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
@@ -42,7 +43,6 @@ void computeNeighboringTetsOfACavity (const std::vector<MTet4*> &cavity,std::vec
}
}
}
-
bool gmshBuildEdgeCavity ( MTet4 *t,
int iLocalEdge,
MVertex **v1,MVertex **v2,
@@ -101,6 +101,9 @@ bool gmshBuildEdgeCavity ( MTet4 *t,
return true;
}
+
+
+// return 1 if it is closed and
typedef struct {
int nbr_triangles ; /* number of different triangles */
int (*triangles)[3] ; /* triangles array */
@@ -260,7 +263,7 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
}
int iBest;
- double best = 0.0;
+ double best = -1.0;
for (int i=0;i<sp.nbr_trianguls;i++){
if(minQuality[i] > best){
best = minQuality[i];
@@ -306,21 +309,6 @@ bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
connectTets ( outside );
-// for (int i=0;i<outside.size();i++){
-// for (int j=0;j<4;j++){
-// MTet4 *neigh = outside[i]->getNeigh(j);
-// // printf("out(%d,%p) neigh(%d) = %p\n",i,outside[i],j,neigh);
-
-// if (neigh && neigh->isDeleted())
-// {
-// bool found = false;
-// for (int k=0;k<cavity.size();k++){
-// if (cavity[k] == neigh) found = true;
-// }
-// printf ("some old deleted tets are still connected to new tets, %d\n",found);
-// }
-// }
-// }
return true;
}
@@ -489,12 +477,95 @@ void gmshBuildVertexCavity_recur ( MTet4 *t,
}
}
+// sliver removal by compound mesh modif
+// postulate : the edge cannot be swopped
+// so we split it, and then collapse the new
+// vertex on another one (of course, not the
+// other one on the unswappable edge)
+// after that crap, the sliver is trashed
+
+bool gmshSliverRemoval ( std::vector<MTet4 *> &newTets,
+ MTet4 *t,
+ const gmshQualityMeasure4Tet &cr){
+ // look if
+ std::vector<MTet4*> cavity;
+ std::vector<MTet4*> outside;
+ std::vector<MVertex*> ring;
+ MVertex *v1,*v2;
+
+ bool isClosed[6];
+ int nbSwappable = 0;
+ int iSwappable;
+ for (int i=0;i<6;i++){
+ isClosed[i] = gmshBuildEdgeCavity ( t,i,&v1,&v2,cavity,outside,ring);
+ if (isClosed[i]){
+ nbSwappable++;
+ iSwappable = i;
+ }
+ }
+
+ if (nbSwappable == 0){
+ // all edges are on model edges or model faces, which means that
+ // nothing can be done
+ return false;
+ }
+ else if (nbSwappable == 1){
+ // classical case, the sliver has 5 edges on the boundary
+ // try to swap first
+ if (gmshEdgeSwap(newTets,t,iSwappable,QMTET_3))return true;
+ // if it does not work, split, smooth and collapse
+ MVertex *v1 = t->tet()->getVertex(edges[iSwappable][0]);
+ MVertex *v2 = t->tet()->getVertex(edges[iSwappable][1]);
+ MVertex *newv = new MVertex (0.5*(v1->x()+v2->x()),
+ 0.5*(v1->y()+v2->y()),
+ 0.5*(v1->z()+v2->z()),t->onWhat());
+ t->onWhat()->mesh_vertices.push_back(newv);
+
+ if (!gmshEdgeSplit(newTets,t,newv,iSwappable,QMTET_ONE))return false;
+ for (int i=0;i<4;i++){
+ if (newTets[newTets.size()-1]->tet()->getVertex(i) == newv)
+ {
+ gmshSmoothVertex(newTets[newTets.size()-1], i,cr);
+ gmshSmoothVertexOptimize (newTets[newTets.size()-1], i,cr);
+ }
+ }
+
+ for (int i=0; i<newTets.size();i++){
+ MTet4 *new_t = newTets[i];
+ if (!(new_t->isDeleted())){
+ for (int j=0;j<6;j++){
+ MVertex *va = new_t->tet()->getVertex(edges[j][0]);
+ MVertex *vb = new_t->tet()->getVertex(edges[j][1]);
+ if (va == newv &&
+ (va != v1 && vb != v1 && va != v2 && vb != v2)){
+ gmshCollapseVertex(newTets,new_t,edges[j][0],edges[j][1],cr);
+ }
+ else if (vb == newv &&
+ (va != v1 && vb != v1 && va != v2 && vb != v2)){
+ gmshCollapseVertex(newTets,new_t,edges[j][1],edges[j][0],cr);
+ }
+ }
+ }
+ }
+
+ return true;
+ }
+ else{
+ for (int i=0;i<4;i++){
+ gmshSmoothVertex(t, i,cr);
+ gmshSmoothVertexOptimize (t, i,cr);
+ }
+ }
+ return false;
+}
bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
MTet4 *t,
int iVertex,
int iTarget,
- const gmshQualityMeasure4Tet &cr){
+ const gmshQualityMeasure4Tet &cr,
+ const gmshLocalMeshModAction action,
+ double *minQual){
if(t->isDeleted())throw;
@@ -512,9 +583,12 @@ bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
std::vector<MTet4*> toDelete;
std::vector<MTet4*> toUpdate;
double volume = 0;
+ double worst = 1.0;
for (int i=0;i<cavity_v.size();i++){
bool found = false;
volume+=fabs(cavity_v[i]->tet()->getVolume());
+ double q = cavity_v[i]->getQuality();
+ worst = std::min(worst,q);
for (int j=0;j<4;j++){
if (cavity_v[i]->tet()->getVertex(j) == tg)found=true;
}
@@ -530,19 +604,30 @@ bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
v->z() = tg->z();
double volume_update=0;
+
+ double worstAfter = 1.0;
+ double newQuals[2000];
+ if (toUpdate.size() >= 2000) throw;
for (int i=0;i<toUpdate.size();i++){
- volume_update+=fabs(toUpdate[i]->tet()->getVolume());
+ double vv;
+ newQuals[i] = qmTet(toUpdate[i]->tet(),cr,&vv);
+ worstAfter = std::min(worstAfter,newQuals[i]);
+ volume_update+=vv;
}
- // printf("%12.5E %12.5E\n",volume,volume_update);
+ // printf("%12.5E %12.5E %12.5E %12.5E %d\n",volume,volume_update,worstAfter,worst,toUpdate.size());
- if (fabs(volume-volume_update) > 1.e-10 * volume)
+ if (fabs(volume-volume_update) > 1.e-10 * volume || worstAfter < worst)
{
v->x() = x;
v->y() = y;
v->z() = z;
return false;
}
+ if (action == GMSH_EVALONLY){
+ *minQual = worstAfter;
+ return true;
+ }
// ok we collapse
computeNeighboringTetsOfACavity (cavity_v,outside);
for (int i=0;i<toUpdate.size();i++){
@@ -552,6 +637,7 @@ bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
toUpdate[i]->tet()->getVertex(3) == v ? tg : toUpdate[i]->tet()->getVertex(3));
MTet4 *t41 = new MTet4 ( tr1 , cr) ;
t41->setOnWhat(cavity_v[0]->onWhat());
+ t41->setQuality(newQuals[i]);
outside.push_back(t41);
newTets.push_back(t41);
}
@@ -562,10 +648,9 @@ bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
return true;
}
-
-bool gmshSmoothVertex ( MTet4 *t,
- int iVertex,
- const gmshQualityMeasure4Tet &cr){
+bool gmshSmoothVertex( MTet4 *t,
+ int iVertex,
+ const gmshQualityMeasure4Tet &cr){
if(t->isDeleted())throw;
if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3)return false;
@@ -611,7 +696,7 @@ bool gmshSmoothVertex ( MTet4 *t,
t->tet()->getVertex(iVertex)->z() = zcg;
double worstAfter = 1.0;
double newQuals[2000];
- if (cavity.size() > 2000) throw;
+ if (cavity.size() >= 2000) throw;
for (int i=0 ; i< cavity.size() ; i++){
double volume;
newQuals[i] = qmTet(cavity[i]->tet(),cr,&volume);
@@ -623,18 +708,128 @@ bool gmshSmoothVertex ( MTet4 *t,
t->tet()->getVertex(iVertex)->x() = x;
t->tet()->getVertex(iVertex)->y() = y;
t->tet()->getVertex(iVertex)->z() = z;
- return false;
+ return false;//gmshSmoothVertexOptimize ( t,iVertex,cr);
}
else{
// restore new quality
for (int i=0 ; i< cavity.size() ; i++){
cavity[i]->setQuality(newQuals[i]);
- }
-
+ }
return true;
}
}
+struct smoothVertexData3D{
+ MVertex *v;
+ std::vector < MTet4 * >ts;
+ double LC;
+};
+
+double smoothing_objective_function_3D(double X, double Y, double Z, MVertex *v, std::vector < MTet4 * >&ts){
+ const double oldX = v->x();
+ const double oldY = v->y();
+ const double oldZ = v->z();
+ v->x() = X;
+ v->y() = Y;
+ v->z() = Z;
+
+ std::vector < MTet4 * >::iterator it = ts.begin();
+ std::vector < MTet4 * >::iterator ite = ts.end();
+ double qMin=1,vol;
+ while(it != ite) {
+ qMin = std::min(qmTet((*it)->tet(),QMTET_2,&vol),qMin);
+ ++it;
+ }
+ v->x() = oldX;
+ v->y() = oldY;
+ v->z() = oldZ;
+ return -qMin;
+}
+
+void deriv_smoothing_objective_function_3D(double X, double Y, double Z,
+ double &F,
+ double &dFdX, double &dFdY, double &dFdZ,
+ void *data){
+ smoothVertexData3D *svd = (smoothVertexData3D*)data;
+ MVertex *v = svd->v;
+ const double LARGE = svd->LC*1.e5;
+ const double SMALL = 1./LARGE;
+ F = smoothing_objective_function_3D(X,Y,Z,v,svd->ts);
+ double F_X = smoothing_objective_function_3D(X+SMALL,Y,Z,v,svd->ts);
+ double F_Y = smoothing_objective_function_3D(X,Y+SMALL,Z,v,svd->ts);
+ double F_Z = smoothing_objective_function_3D(X,Y,Z+SMALL,v,svd->ts);
+ dFdX = (F_X-F)*LARGE;
+ dFdY = (F_Y-F)*LARGE;
+ dFdZ = (F_Z-F)*LARGE;
+}
+
+double smooth_obj_3D(double X, double Y, double Z, void *data){
+ smoothVertexData3D *svd = (smoothVertexData3D*)data;
+ return smoothing_objective_function_3D(X,Y,Z,svd->v,svd->ts);
+}
+
+
+bool gmshSmoothVertexOptimize ( MTet4 *t,
+ int iVertex,
+ const gmshQualityMeasure4Tet &cr){
+
+ if(t->tet()->getVertex(iVertex)->onWhat()->dim() < 3)return false;
+
+ smoothVertexData3D vd;
+ vd.ts.push_back(t);
+ vd.v = t->tet()->getVertex(iVertex);
+ vd.LC = 1.0; // WRONG
+ gmshBuildVertexCavity_recur (t,t->tet()->getVertex(iVertex),vd.ts);
+
+ double xopti=vd.v->x();
+ double yopti=vd.v->y();
+ double zopti=vd.v->z();
+
+ double val;
+ minimize_3 ( smooth_obj_3D, deriv_smoothing_objective_function_3D, &vd, 4, xopti,yopti,zopti,val);
+
+ double vTot=0;
+
+ for (int i=0 ; i< vd.ts.size() ; i++){
+ double volume = fabs(vd.ts[i]->tet()->getVolume());
+ vTot += volume;
+ }
+
+ double volumeAfter = 0.0;
+
+ double x = t->tet()->getVertex(iVertex)->x();
+ double y = t->tet()->getVertex(iVertex)->y();
+ double z = t->tet()->getVertex(iVertex)->z();
+
+ t->tet()->getVertex(iVertex)->x() = xopti;
+ t->tet()->getVertex(iVertex)->y() = yopti;
+ t->tet()->getVertex(iVertex)->z() = zopti;
+
+ double newQuals[2000];
+ if (vd.ts.size() >= 2000) throw;
+ for (int i=0 ; i< vd.ts.size() ; i++){
+ double volume;
+ newQuals[i] = qmTet(vd.ts[i]->tet(),cr,&volume);
+ volumeAfter += volume;
+ }
+
+ if (fabs(volumeAfter-vTot) > 1.e-10 * vTot){
+ t->tet()->getVertex(iVertex)->x() = x;
+ t->tet()->getVertex(iVertex)->y() = y;
+ t->tet()->getVertex(iVertex)->z() = z;
+ return false;
+ }
+ else{
+ // restore new quality
+ for (int i=0 ; i< vd.ts.size() ; i++){
+ vd.ts[i]->setQuality(newQuals[i]);
+ }
+ return true;
+ }
+}
+
+
+
// Edge split sets ...
// Here, we only build a list of tets that are a subdivision
diff --git a/Mesh/meshGRegionLocalMeshMod.h b/Mesh/meshGRegionLocalMeshMod.h
index 1342d6038132bc052b0b920a7accd6e0e962c02f..af77af2b128794aee34efae4f56ed11c106173b5 100644
--- a/Mesh/meshGRegionLocalMeshMod.h
+++ b/Mesh/meshGRegionLocalMeshMod.h
@@ -26,6 +26,9 @@
// operators only apply to the "bulk" of the
// mesh and cannot be applied to boudnaries.
// I'm working on it
+
+enum gmshLocalMeshModAction {GMSH_DOIT,GMSH_EVALONLY};
+
bool gmshEdgeSwap (std::vector<MTet4 *> &newTets,
MTet4 *tet,
int iLocalEdge,
@@ -40,11 +43,17 @@ bool gmshSmoothVertex ( MTet4 *t,
int iLocalVertex,
const gmshQualityMeasure4Tet &cr);
+bool gmshSmoothVertexOptimize ( MTet4 *t,
+ int iVertex,
+ const gmshQualityMeasure4Tet &cr);
+
bool gmshCollapseVertex ( std::vector<MTet4 *> &newTets,
MTet4 *t,
int iVertex,
int iTarget,
- const gmshQualityMeasure4Tet &cr);
+ const gmshQualityMeasure4Tet &cr,
+ const gmshLocalMeshModAction = GMSH_DOIT,
+ double *result = 0);
bool gmshEdgeSplit (std::vector<MTet4 *> &newTets,
MTet4 *tet,
@@ -52,6 +61,9 @@ bool gmshEdgeSplit (std::vector<MTet4 *> &newTets,
int iLocalEdge,
const gmshQualityMeasure4Tet &cr);
+bool gmshSliverRemoval ( std::vector<MTet4 *> &newTets,
+ MTet4 *t,
+ const gmshQualityMeasure4Tet &cr);
#endif
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 10530237031ba4c6f540e7577c1f68fda0ae25ba..5f6e2a9e32422e71c8898f96988e60a6a2573e14 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -84,6 +84,8 @@ double qmTet ( const double &x1, const double &y1, const double &z1,
double quality;
switch(cr)
{
+ case QMTET_ONE:
+ return 1.0;
case QMTET_3:
{
double mat[3][3];
diff --git a/Mesh/qualityMeasures.h b/Mesh/qualityMeasures.h
index 6b4a974e03365aeccb159aecba3267bc5785bdf4..96d9078bd012fbd02e35a9052045dabc6f5dbc18 100644
--- a/Mesh/qualityMeasures.h
+++ b/Mesh/qualityMeasures.h
@@ -7,7 +7,7 @@ class MVertex;
class MTriangle;
class MTetrahedron;
enum gmshQualityMeasure4Triangle {QMTRI_RHO};
-enum gmshQualityMeasure4Tet {QMTET_1,QMTET_2,QMTET_3};
+enum gmshQualityMeasure4Tet {QMTET_1,QMTET_2,QMTET_3,QMTET_ONE};
double qmTriangle ( MTriangle *f, const gmshQualityMeasure4Triangle &cr);
double qmTriangle ( BDS_Face *f, const gmshQualityMeasure4Triangle &cr);
diff --git a/Numeric/Makefile b/Numeric/Makefile
index 58705f7fa6f6018b333ad8746308e21dc1661fd0..d688c90139f3ac8de26b98e99e14df6d4addd57d 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.37 2007-11-12 10:49:16 geuzaine Exp $
+# $Id: Makefile,v 1.38 2008-01-14 21:29:14 remacle Exp $
#
# Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
#
@@ -29,6 +29,7 @@ SRC = Numeric.cpp\
EigSolve.cpp\
predicates.cpp\
gsl_newt.cpp\
+ gsl_min.cpp\
gsl_brent.cpp
OBJ = ${SRC:.cpp=.o}
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index b387f416956192c5a8b9dd23ab0b374db2ef185b..d87c3337f6bc4afb21c678616d132d4ff2c5c7ca 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -1,4 +1,4 @@
-// $Id: Numeric.cpp,v 1.36 2007-11-11 19:53:57 remacle Exp $
+// $Id: Numeric.cpp,v 1.37 2008-01-14 21:29:14 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -236,6 +236,16 @@ double det2x2(double mat[2][2])
return mat[0][0] *mat[1][1] -mat[1][0] *mat[0][1];
}
+double det2x3(double mat[2][3])
+{
+ double v1[3] = {mat[0][0],mat[0][1],mat[0][2]};
+ double v2[3] = {mat[1][0],mat[1][1],mat[1][2]};
+ double n[3];
+
+ prodve (v1,v2,n);
+ return norm3(n);
+}
+
double inv2x2(double mat[2][2], double inv[2][2])
{
const double det = det2x2 ( mat );
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index da6d16ef5889fe7fd7978a1caae26c45b52b412c..326a1a6b9e5088fd37c97f1671c501e14a773b07 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -90,6 +90,7 @@ int sys2x2(double mat[2][2], double b[2], double res[2]);
int sys3x3(double mat[3][3], double b[3], double res[3], double *det);
int sys3x3_with_tol(double mat[3][3], double b[3], double res[3], double *det);
double det2x2(double mat[2][2]);
+double det2x3(double mat[2][3]);
double det3x3(double mat[3][3]);
double trace3x3(double mat[3][3]);
double trace2 (double mat[3][3]);
@@ -121,6 +122,19 @@ void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb,
double *fc, double (*func)(double));
void newt(double x[], int n, int *check,
void (*vecfunc)(int, double [], double []));
+void minimize_2 ( double (*f) (double, double, void *data),
+ void (*df) (double, double, double &, double &, double &, void *data) ,
+ void *data,int niter,
+ double &U, double &V, double &res);
+void minimize_3 ( double (*f) (double, double, double, void *data),
+ void (*df) (double, double, double , double &, double &, double &, double &, void *data) ,
+ void *data,int niter,
+ double &U, double &V, double &W, double &res);
+void minimize_N (int N,
+ double (*f) (double*, void *data),
+ void (*df) (double*, double*, double &, void *data) ,
+ void *data,int niter,
+ double *, double &res);
/* Robust geometrical predicates */
namespace gmsh{
diff --git a/benchmarks/2d/coin.geo b/benchmarks/2d/coin.geo
index 0f596c122bd81a1c51469df227fdb5c1dfbba785..8ec850f234bda98f1436aafcc369c56bd40dcd04 100644
--- a/benchmarks/2d/coin.geo
+++ b/benchmarks/2d/coin.geo
@@ -4,14 +4,14 @@ a = 1.0;
b = 1.0;
h = 1.0;
-refine = 0.0001;
+refine = 0.0000001;
Point(1) = {0, 0, 0, lc*refine};
Point(2) = {a, 0, 0, lc} ;
Point(3) = {0, b, 0, lc} ;
Point(4) = {a, -h, 0, lc} ;
Point(5) = {-h, b, 0, lc} ;
-Point(6) = {-h, -h, 0, lc} ;
+Point(6) = {-h, -h, 0, lc*refine} ;
Line(1) = {1,3} ;
Line(2) = {3,5} ;
diff --git a/benchmarks/2d/function.geo b/benchmarks/2d/function.geo
index 0d1f9f4cbdd0122fac6dcf1f7d39435985db5f2a..896c924d9d45b2b16622afa75581bcf943fbdc6a 100644
--- a/benchmarks/2d/function.geo
+++ b/benchmarks/2d/function.geo
@@ -5,15 +5,15 @@ theloop = 0;
Function myCircle
p1 = newp;
- Point (p1) = {x,y,0,0.2};
+ Point (p1) = {x,y,0,0.4};
p2 = newp;
- Point (p2) = {r+x,y,0,0.2};
+ Point (p2) = {r+x,y,0,0.4};
p3 = newp;
- Point (p3) = {x,r+y,0,0.2};
+ Point (p3) = {x,r+y,0,0.4};
p4 = newp;
- Point (p4) = {-r+x,y,0,0.2};
+ Point (p4) = {-r+x,y,0,0.4};
p5 = newp;
- Point (p5) = {x,-r+y,0,0.2};
+ Point (p5) = {x,-r+y,0,0.4};
c1 = newreg;
Circle (c1) = {p2,p1,p3};
c2 = newreg;
@@ -51,7 +51,5 @@ Call myCircle;
loop5 = theloop;
Plane Surface(newreg) = {loop5,loop4,loop3,loop2,loop1};
-Line(10000) = {6,11};
-Attractor Line {10000} = {0.5,0.05,5,5,100};
diff --git a/benchmarks/2d/wing-splines.geo b/benchmarks/2d/wing-splines.geo
index 1f2c8da7d420e81f45b6980ffed8c93c08e2c8ea..83a359af18a664a0d4a263ffaa1f309f4a3bad62 100644
--- a/benchmarks/2d/wing-splines.geo
+++ b/benchmarks/2d/wing-splines.geo
@@ -624,7 +624,7 @@ Circle(409) = {4354,4351,4353};
Circle(410) = {4353,4351,4352};
// use an embedded curve
-Line {407,408,409,410} In Surface {406};
+//Line {407,408,409,410} In Surface {406};
// or an attractor
//Attractor Line {1 ... 12} = {1,lc_wing,lc_box,1000,3};