From 147ed21a917aba459afaa7e048cf97eadc2793e4 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 16 Mar 2007 10:03:40 +0000
Subject: [PATCH] *** empty log message ***
---
Geo/MElement.cpp | 10 ++---
Geo/OCCFace.cpp | 6 +--
Mesh/Attractors.cpp | 18 ++++-----
Mesh/HighOrder.cpp | 98 +++++++++++++++++++++++++--------------------
Mesh/HighOrder.h | 1 +
Mesh/meshGFace.cpp | 9 +++--
Parser/OpenFile.cpp | 5 ++-
7 files changed, 81 insertions(+), 66 deletions(-)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 650a2fbd0f..5d05978a34 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.32 2007-03-14 15:47:49 remacle Exp $
+// $Id: MElement.cpp,v 1.33 2007-03-16 10:03:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -548,8 +548,8 @@ void GradGeomShapeFunctionP4 (double u, double v, double grads[12][2])
grads[i][0] = 0;
grads[i][1] = 0;
for (int j=0;j<12;j++){
- if (P4[j][0] > 0)grads[i][0] += coef3[i][j] * pow(u,P4[j][0] - 1 ) * pow(v,P4[j][1] ) ;
- if (P4[j][1] > 0)grads[i][1] += coef3[i][j] * pow(u,P4[j][0] ) * pow(v,P4[j][1] -1 ) ;
+ 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 ) ;
}
}
}
@@ -560,8 +560,8 @@ 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] += coef3[i][j] * pow(u,P5[j][0] - 1 ) * pow(v,P5[j][1] ) ;
- if (P5[j][1] > 0)grads[i][1] += coef3[i][j] * pow(u,P5[j][0] ) * pow(v,P5[j][1] -1 ) ;
+ 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 ) ;
}
}
}
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 28a95892cd..06d41e59b5 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.21 2007-02-27 17:15:47 remacle Exp $
+// $Id: OCCFace.cpp,v 1.22 2007-03-16 10:03:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -41,7 +41,7 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
TopExp_Explorer exp0, exp01, exp1, exp2, exp3;
for(exp2.Init(s, TopAbs_WIRE); exp2.More(); exp2.Next()){
TopoDS_Shape wire = exp2.Current();
- Msg(INFO,"OCC Face %d - New Wire",num);
+ Msg(DEBUG2,"OCC Face %d - New Wire",num);
std::list<GEdge*> l_wire;
for(exp3.Init(wire, TopAbs_EDGE); exp3.More(); exp3.Next()){
TopoDS_Edge edge = TopoDS::Edge(exp3.Current());
@@ -70,7 +70,7 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num, TopTools_IndexedMapOfShape
_periodic[1] = surface.IsVPeriodic();
ShapeAnalysis::GetFaceUVBounds(_s, umin, umax, vmin, vmax);
- Msg(INFO, "OCC Face %d with %d edges bounds (%g,%g)(%g,%g)", num, l_edges.size(),umin,umax,vmin,vmax);
+ Msg(DEBUG2, "OCC Face %d with %d edges bounds (%g,%g)(%g,%g)", num, l_edges.size(),umin,umax,vmin,vmax);
// we do that for the projections to converge on the
// borders of the surface
const double du = umax-umin;
diff --git a/Mesh/Attractors.cpp b/Mesh/Attractors.cpp
index c77c7dee14..ee365b7ebe 100644
--- a/Mesh/Attractors.cpp
+++ b/Mesh/Attractors.cpp
@@ -1,4 +1,4 @@
-// $Id: Attractors.cpp,v 1.3 2007-02-27 17:15:47 remacle Exp $
+// $Id: Attractors.cpp,v 1.4 2007-03-16 10:03:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -69,21 +69,21 @@ void Attractor::buildFastSearchStructures()
delete kdtree;
}
int totpoints = 0;
- for(std::list <Attractor *>::iterator it = allAttractors.begin();
- it != allAttractors.end(); ++it)
- totpoints += (*it)->attractorPoints.size();
+// for(std::list <Attractor *>::iterator it = allAttractors.begin();
+// it != allAttractors.end(); ++it)
+ totpoints += attractorPoints.size();
if (totpoints)
zeronodes = annAllocPts(totpoints, 4);
int k = 0;
- for(std::list <Attractor *>::iterator it = allAttractors.begin();
- it != allAttractors.end(); ++it){
- for (std::list <SPoint3>::iterator it2 = (*it)->attractorPoints.begin();
- it2 != (*it)->attractorPoints.end(); ++it2){
+// for(std::list <Attractor *>::iterator it = allAttractors.begin();
+// it != allAttractors.end(); ++it){
+ for (std::list <SPoint3>::iterator it2 = attractorPoints.begin();
+ it2 != attractorPoints.end(); ++it2){
zeronodes[k][0]=it2->x();
zeronodes[k][1]=it2->y();
zeronodes[k++][2]=it2->z();
}
- }
+// }
kdtree = new ANNkd_tree(zeronodes, totpoints, 3);
#endif
}
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 4d76f08a20..f215c3a021 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.9 2007-03-14 15:47:49 remacle Exp $
+// $Id: HighOrder.cpp,v 1.10 2007-03-16 10:03:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -631,12 +631,19 @@ void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
double u[7] = {0,1,0,.5,0,.5,.3333};
double v[7] = {0,0,1,.5,.5,0,.3333};
double j[2][2];
- for (int i=0;i<7;i++)
+
+ int n = 3;
+
+ for (int i=0;i<n;i++)
{
- t->jac ( u[i], v[i] , j );
- double det = det2x2 ( j );
- minJ = std::min ( det , minJ );
- maxJ = std::max ( det , maxJ );
+ for (int k=0;k<n-i;k++)
+ {
+ t->jac ( (double)i/(n-1), (double)k/(n-1) , j );
+ double det = det2x2 ( j );
+ // printf("det = %12.5E\n",det);
+ minJ = std::min ( det , minJ );
+ maxJ = std::max ( det , maxJ );
+ }
}
}
@@ -742,11 +749,47 @@ void smoothInternalEdges ( GFace *gf , edgeContainer &edgeVertices)
}
}
- }
-
-
+ }
}
+void checkHighOrderTriangles ( GModel *m )
+{
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
+ {
+ bool twod = true;
+// for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
+// 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)
+ {
+ 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(GERROR,"Triangle %d %d %d has negative jacobians in GFace %d",t->getVertex(0)->getNum(),
+ t->getVertex(1)->getNum(),t->getVertex(2)->getNum(),(*it)->tag());
+ }
+ }
+ if (minJ != 1.e22)
+ {
+ if (minJ*maxJ < 0)
+ Msg(GERROR,"There exists triangles with negative jacobians in GFace %d",(*it)->tag());
+ else
+ Msg(INFO,"No negative jacobians detected on model face %d range = (%g,%g)",(*it)->tag(),minJ,maxJ);
+ }
+ }
+ }
+}
+
void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
{
@@ -789,43 +832,10 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
if (CTX.mesh.smooth_internal_edges)
for (int i=0;i<10;i++)
smoothInternalEdges ( *it , edgeVertices);
+ }
- bool twod = true;
- for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
- 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;
- }
-
+ checkHighOrderTriangles ( m );
- if (twod && order > 3)
- {
- Msg(INFO,"Validity check of higher order meshes is still bugged for orders > 3");
- }
- else if (twod)
- {
- 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];
- getMinMaxJac (t,minJloc,maxJloc);
- minJ = std::min(minJ,minJloc);
- maxJ = std::max(maxJ,maxJloc);
- if (minJloc*maxJloc < 0)
- Msg(GERROR,"Triangle %d %d %d has negative jacobians in GFace %d",t->getVertex(0)->getNum(),
- t->getVertex(1)->getNum(),t->getVertex(2)->getNum(),(*it)->tag());
- }
-
- if (minJ*maxJ < 0)
- Msg(GERROR,"There exists triangles with negative jacobians in GFace %d",(*it)->tag());
- else
- Msg(INFO,"Jacobian range of face %d is %g %g",(*it)->tag(),minJ,maxJ);
-
- }
- }
-
double t2 = Cpu();
Msg(INFO, "Mesh second order complete (%g s)", t2 - t1);
Msg(STATUS1, "Mesh");
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index 84815275f8..d81772dc98 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -24,5 +24,6 @@
void SetOrder1(GModel *m);
void SetOrderN(GModel *m, int order, bool linear=true, bool incomplete=false);
+void checkHighOrderTriangles ( GModel *m );
#endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 8bf9dcfea3..0fa7c0bb0f 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.67 2007-03-14 11:44:18 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.68 2007-03-16 10:03:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -873,8 +873,9 @@ bool gmsh2DMeshGenerator ( GFace *gf , bool debug = true)
// if (CTX.mesh.algo2d == ALGO_2D_MESHADAPT)
{
RefineMesh (gf,*m,10);
- OptimizeMesh(gf, *m, 3);
+ OptimizeMesh(gf, *m, 2);
RefineMesh (gf,*m,-10);
+ OptimizeMesh(gf, *m, 2);
if (gf->meshAttributes.recombine)
{
m->recombineIntoQuads (gf->meshAttributes.recombineAngle,gf);
@@ -1403,9 +1404,9 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf , bool debug = true)
// if (CTX.mesh.algo2d == ALGO_2D_MESHADAPT)
{
RefineMesh (gf,*m,10);
- OptimizeMesh(gf, *m, 5);
+ OptimizeMesh(gf, *m, 2);
RefineMesh (gf,*m,-10);
- OptimizeMesh(gf, *m, 1);
+ OptimizeMesh(gf, *m, 2);
if (gf->meshAttributes.recombine)
{
diff --git a/Parser/OpenFile.cpp b/Parser/OpenFile.cpp
index 2680a7b672..7474acff9d 100644
--- a/Parser/OpenFile.cpp
+++ b/Parser/OpenFile.cpp
@@ -1,4 +1,4 @@
-// $Id: OpenFile.cpp,v 1.144 2007-03-11 20:19:05 geuzaine Exp $
+// $Id: OpenFile.cpp,v 1.145 2007-03-16 10:03:40 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -34,6 +34,7 @@
#include "Views.h"
#include "ReadImg.h"
#include "OS.h"
+#include "HighOrder.h"
#if defined(HAVE_FLTK)
#include "GmshUI.h"
@@ -396,6 +397,8 @@ int MergeFile(char *name, int warn_if_missing)
CTX.geom.draw = 1;
CTX.mesh.changed = ENT_ALL;
+ checkHighOrderTriangles ( GMODEL );
+
Msg(STATUS2, "Read '%s'", name);
return status;
}
--
GitLab