From 510c88ead55e65374dec01937dc90c03f6c45339 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 8 Jul 2008 12:44:34 +0000
Subject: [PATCH] mesh drawing debugged INRIA input corrected and enhanced
---
Geo/GFace.cpp | 18 +++++++++++--
Geo/GModelIO_Mesh.cpp | 56 +++++++++++++++++++++++++++++++++++++----
Geo/MElement.cpp | 26 ++++++++++---------
Geo/OCCFace.cpp | 4 +--
Geo/gmshFace.cpp | 21 +---------------
Geo/gmshFace.h | 1 -
Mesh/HighOrder.cpp | 4 +--
Numeric/FunctionSpace.h | 1 -
8 files changed, 86 insertions(+), 45 deletions(-)
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 453172c067..a450c243c1 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.67 2008-07-04 12:03:50 geuzaine Exp $
+// $Id: GFace.cpp,v 1.68 2008-07-08 12:44:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -730,7 +730,21 @@ bool GFace::buildSTLTriangulation()
// by default we assume that straight lines are geodesics
SPoint2 GFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
{
- return pt1 + (pt2 - pt1) * t;
+ if(CTX.mesh.second_order_experimental){
+ // FIXME: this is buggy -- remove the CTX option once we do it in
+ // a robust manner
+ GPoint gp1 = point(pt1.x(), pt1.y());
+ GPoint gp2 = point(pt2.x(), pt2.y());
+ SPoint2 guess = pt1 + (pt2 - pt1) * t;
+ GPoint gp = closestPoint(SPoint3(gp1.x() + t * (gp2.x() - gp1.x()),
+ gp1.y() + t * (gp2.y() - gp1.y()),
+ gp1.z() + t * (gp2.z() - gp1.z())),
+ (double*)guess);
+ return SPoint2(gp.u(), gp.v());
+ }
+ else{
+ return pt1 + (pt2 - pt1) * t;
+ }
}
// length of a curve drawn on a surface
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index cf7d039555..e13d701e10 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.57 2008-07-04 16:58:48 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.58 2008-07-08 12:44:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -1396,18 +1396,19 @@ int GModel::readMESH(const std::string &name)
char str[256];
int format;
sscanf(buffer, "%s %d", str, &format);
-
if(format != 1){
Msg::Error("Medit mesh import only available for ASCII files");
return 0;
}
std::vector<MVertex*> vertexVector;
- std::map<int, std::vector<MElement*> > elements[3];
+ std::map<int, std::vector<MElement*> > elements[4];
+ std::vector<MVertex*> corners,ridges;
while(!feof(fp)) {
- if(!fgets(buffer, sizeof(buffer), fp)) break;
- if(buffer[0] != '#'){ // skip comments
+ if(!fgets(buffer, 256, fp)) break;
+ if(buffer[0] != '#'){ // skip comments and empty lines
+ str[0]='\0';
sscanf(buffer, "%s", str);
if(!strcmp(str, "Dimension")){
if(!fgets(buffer, sizeof(buffer), fp)) break;
@@ -1426,6 +1427,21 @@ int GModel::readMESH(const std::string &name)
vertexVector[i] = new MVertex(x, y, z);
}
}
+ else if(!strcmp(str, "Edges")){
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int nbe;
+ sscanf(buffer, "%d", &nbe);
+ Msg::Info("%d edges", nbe);
+ for(int i = 0; i < nbe; i++) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[2], cl;
+ sscanf(buffer, "%d %d %d", &n[0], &n[1], &cl);
+ for(int j = 0; j < 2; j++) n[j]--;
+ std::vector<MVertex*> vertices;
+ if(!getVertices(2, n, vertexVector, vertices)) return 0;
+ elements[3][cl].push_back(new MLine(vertices));
+ }
+ }
else if(!strcmp(str, "Triangles")){
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
@@ -1441,6 +1457,36 @@ int GModel::readMESH(const std::string &name)
elements[0][cl].push_back(new MTriangle(vertices));
}
}
+ else if(!strcmp(str, "Corners")){
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int nbe;
+ sscanf(buffer, "%d", &nbe);
+ Msg::Info("%d corners", nbe);
+ for(int i = 0; i < nbe; i++) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[1];
+ sscanf(buffer, "%d",&n[0]);
+ for(int j = 0; j < 1; j++) n[j]--;
+ // std::vector<MVertex*> vertices;
+ // if(!getVertices(1, n, vertexVector, vertices)) return 0;
+ // corners.push_back(vertices[0]);
+ }
+ }
+ else if(!strcmp(str, "Ridges")){
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int nbe;
+ sscanf(buffer, "%d", &nbe);
+ Msg::Info("%d ridges", nbe);
+ for(int i = 0; i < nbe; i++) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[1];
+ sscanf(buffer, "%d",&n[0]);
+ for(int j = 0; j < 1; j++) n[j]--;
+ // std::vector<MVertex*> vertices;
+ // if(!getVertices(1, n, vertexVector, vertices)) return 0;
+ // ridges.push_back(vertices[0]);
+ }
+ }
else if(!strcmp(str, "Quadrilaterals")) {
if(!fgets(buffer, sizeof(buffer), fp)) break;
int nbe;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 4e85bfe43f..8b7b456782 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.77 2008-07-03 17:06:01 geuzaine Exp $
+// $Id: MElement.cpp,v 1.78 2008-07-08 12:44:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -690,27 +690,29 @@ void MTriangle::pnt(int ord, MVertex *vs[], double uu, double vv, double ww, SPo
double sf[256];
int nf = getNumFaceVertices();
+ // printf("%d %d\n",nf,ord);
if (!nf){
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, ww,sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, ww,sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv, ww,sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv, ww,sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv, ww,sf); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv,sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv,sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_9).f(uu, vv,sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_12).f(uu, vv,sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_15I).f(uu, vv,sf); break;
default: Msg::Error("Order %d triangle pnt not implemented", ord); break;
}
}
else{
switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv, ww,sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv, ww,sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv, ww,sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv, ww,sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv, ww,sf); break;
+ case 1: gmshFunctionSpaces::find(MSH_TRI_3).f(uu, vv,sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TRI_6).f(uu, vv,sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TRI_10).f(uu, vv,sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TRI_15).f(uu, vv,sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TRI_21).f(uu, vv,sf); break;
default: Msg::Error("Order %d triangle pnt not implemented", ord); break;
}
}
-
+ // printf("coucou\n");
+
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();
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index bab19e83c1..7e22199e88 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -1,4 +1,4 @@
-// $Id: OCCFace.cpp,v 1.43 2008-07-03 17:06:02 geuzaine Exp $
+// $Id: OCCFace.cpp,v 1.44 2008-07-08 12:44:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -157,7 +157,7 @@ GPoint OCCFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) c
double pp[2] = {initialGuess[0],initialGuess[1]};
proj.LowerDistanceParameters(pp[0], pp[1]);
- Msg::Info("projection lower distance parameters %g %g",pp[0],pp[1]);
+ // Msg::Info("projection lower distance parameters %g %g",pp[0],pp[1]);
if((pp[0] < umin || umax < pp[0]) || (pp[1]<vmin || vmax<pp[1])){
Msg::Error("Point projection is out of face bounds");
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 7ec6a112a6..3a0f8edd24 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.65 2008-07-03 17:06:02 geuzaine Exp $
+// $Id: gmshFace.cpp,v 1.66 2008-07-08 12:44:33 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -248,25 +248,6 @@ GEntity::GeomType gmshFace::geomType() const
}
}
-// by default we assume that straight lines are geodesics
-SPoint2 gmshFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
-{
- if(isSphere && CTX.mesh.second_order_experimental){
- // FIXME: this is buggy -- remove the CTX option once we do it in
- // a robust manner
- GPoint gp1 = point(pt1.x(), pt1.y());
- GPoint gp2 = point(pt2.x(), pt2.y());
- SPoint2 guess = pt1 + (pt2 - pt1) * t;
- GPoint gp = closestPoint(SPoint3(gp1.x() + t * (gp2.x() - gp1.x()),
- gp1.y() + t * (gp2.y() - gp1.y()),
- gp1.z() + t * (gp2.z() - gp1.z())),
- (double*)guess);
- return SPoint2(gp.u(), gp.v());
- }
- else{
- return pt1 + (pt2 - pt1) * t;
- }
-}
bool gmshFace::containsPoint(const SPoint3 &pt) const
{
diff --git a/Geo/gmshFace.h b/Geo/gmshFace.h
index e9900b4f1c..c18a3b67da 100644
--- a/Geo/gmshFace.h
+++ b/Geo/gmshFace.h
@@ -34,7 +34,6 @@ class gmshFace : public GFace {
virtual ~gmshFace(){}
Range<double> parBounds(int i) const;
void setModelEdges(std::list<GEdge*>&);
- virtual SPoint2 geodesic(const SPoint2 &pt1, const SPoint2 &pt2, double t);
virtual GPoint point(double par1, double par2) const;
virtual GPoint closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const;
virtual bool containsPoint(const SPoint3 &pt) const;
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index a67d10c4d1..1363055c2d 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.32 2008-07-03 17:06:03 geuzaine Exp $
+// $Id: HighOrder.cpp,v 1.33 2008-07-08 12:44:34 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -1358,7 +1358,7 @@ void printJacobians(GModel *m, const char *nm)
{
return;
- const int n = 5;
+ const int n = 15;
double D[n][n];
double X[n][n];
double Y[n][n];
diff --git a/Numeric/FunctionSpace.h b/Numeric/FunctionSpace.h
index 267347278d..b349898444 100644
--- a/Numeric/FunctionSpace.h
+++ b/Numeric/FunctionSpace.h
@@ -36,7 +36,6 @@ struct gmshFunctionSpace
p[j][1] = pow(vv,monomials(j, 1));
}
}
-
inline void f(double u, double v, double w, double *sf) const
{
for(int i = 0; i < coefficients.size1(); i++){
--
GitLab