From 1ced34ea05c206b5e504957813f43876dbc5d2be Mon Sep 17 00:00:00 2001
From: Koen Hillewaert <koen.hillewaert@cenaero.be>
Date: Mon, 2 Feb 2009 12:36:13 +0000
Subject: [PATCH] a posteriori snapping of higher-order points on geometry
---
Common/OpenFile.cpp | 8 +++
Geo/GEdge.cpp | 111 ++++++++++++++++++++++++++++++++++++++++++
Geo/GEdge.h | 30 +++++++++++-
Geo/GModelIO_Mesh.cpp | 12 +++++
Geo/MVertex.cpp | 14 ++++--
Geo/MVertex.h | 3 +-
Geo/discreteEdge.cpp | 4 ++
Geo/discreteEdge.h | 1 +
Mesh/HighOrder.cpp | 41 ++++++++++------
9 files changed, 205 insertions(+), 19 deletions(-)
diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index c9a7c4b5e3..2e14e19279 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -31,6 +31,10 @@
#include "Draw.h"
#endif
+// KH: modify grid order after reading
+extern Context_T CTX;
+extern void SetOrderN(GModel *m, int order, bool linear, bool incomplete);
+// end KH
#define SQU(a) ((a)*(a))
static void FinishUpBoundingBox()
@@ -346,6 +350,10 @@ int MergeFile(std::string fileName, bool warnIfMissing)
#if !defined(HAVE_NO_POST)
if(status > 1) status = PView::readMSH(fileName);
#endif
+ // KH - modify mesh order after reading
+ if (CTX.mesh.order > 1) SetOrderN(GModel::current(),CTX.mesh.order,false,false);
+ // end - KH
+
}
#if !defined(HAVE_NO_POST)
else if(!strncmp(header, "$PostFormat", 11) ||
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index a59b67bb14..44b3a04127 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -5,6 +5,8 @@
#include <sstream>
#include <algorithm>
+#include <limits>
+
#include "GmshConfig.h"
#include "GmshDefines.h"
#include "GmshMessage.h"
@@ -15,8 +17,11 @@
#if !defined(HAVE_GMSH_EMBEDDED)
#include "GaussLegendre1D.h"
+#include "Context.h"
#endif
+extern Context_T CTX;
+
GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
: GEntity(model, tag), _tooSmall(false), v0(_v0), v1(_v1)
{
@@ -167,6 +172,7 @@ double GEdge::curvature(double par) const
const double one_over_D = 1. / D;
SVector3 d = one_over_D * (n2 - n1);
return norm(d);
+
}
double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints)
@@ -188,3 +194,108 @@ double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints)
#endif
}
+GPoint GEdge::closestPoint(const SPoint3 & q,double& t) const {
+
+ double tolerance = 1.e-12;
+ double dist = 1.;
+
+ Range<double> interval = parBounds(0);
+
+ double tMin = std::min(interval.high(),interval.low());
+ double tMax = std::max(interval.high(),interval.low());
+ double relax = 1.;
+ double dt,dt0,t0;
+ int nb = 10;
+
+ t = (tMin + tMax) * 0.5;
+
+ while (relax > 0.1) {
+
+ int i = 0;
+
+ t = 0.5 * (tMin + tMax);
+ dt0 = tMax - tMin;
+ dt = dt0;
+
+ while (dt > tolerance * dt0 && i++ < nb) {
+ t0 = t;
+ dt0 = dt;
+ SVector3 dp = q - position(t);
+ SVector3 derP = firstDer(t);
+ double b = dot(derP,derP);
+ double c = dot(derP,dp);
+ dt = c / b;
+ t = std::max(tMin,std::min(tMax,t0+relax * dt));
+ dt = fabs(t - t0);
+ }
+ if (i > nb) relax *= 0.5;
+ }
+ return point(t);
+}
+
+
+#include <iostream>
+
+
+double GEdge::parFromPoint(const SVector3& Q) const {
+ double t;
+ bool success = XYZToU(Q,t);
+ return t;
+}
+
+bool GEdge::XYZToU(const SVector3& Q,double &u, const double relax) const
+{
+
+
+ const double Precision = 1.e-8;
+ const int MaxIter = 25;
+ const int NumInitGuess = 11;
+
+ double err, err2;
+ int iter;
+
+ Range<double> uu = parBounds(0);
+ double uMin = uu.low();
+ double uMax = uu.high();
+
+ SVector3 P;
+
+ double init[NumInitGuess];
+
+ for (int i=0;i<NumInitGuess;i++) init[i] = uMin + (uMax-uMin)/(NumInitGuess-1)*i;
+
+ for(int i = 0; i < NumInitGuess; i++){
+
+ u = init[i];
+ double uNew = u;
+ err = 1.0;
+ iter = 1;
+
+ SVector3 dPQ = P - Q;
+ err2 = dPQ.norm();
+
+ if (err2 < 1.e-8 * CTX.lc) return true;
+
+ while(iter++ < MaxIter && err2 > 1e-8 * CTX.lc) {
+ SVector3 der = firstDer(u);
+ uNew = u - relax * dot(dPQ,der) / dot(der,der);
+ uNew = std::min(uMax,std::max(uMin,uNew));
+ P = position(uNew);
+
+ dPQ = P - Q;
+ err2 = dPQ.norm();
+ err = fabs(uNew - u);
+ u = uNew;
+ }
+
+ if (err2 < 1e-8 * CTX.lc) return true;
+ }
+
+ if(relax > 1.e-2) {
+ Msg::Info("point %g %g %g on edge %d : Relaxation factor = %g", Q.x(), Q.y(), Q.z(), 0.75 * relax);
+ return XYZToU(Q, u, 0.75 * relax);
+ }
+
+ Msg::Error("Could not converge reparametrisation of point (%e,%e,%e) on edge %d",Q.x(),Q.y(),Q.z(),tag());
+ return false;
+}
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index 0ac4cfc769..444847c557 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -16,6 +16,8 @@ class MElement;
class MLine;
class ExtrudeParams;
+#include <set>
+
// A model edge.
class GEdge : public GEntity {
private:
@@ -26,6 +28,8 @@ class GEdge : public GEntity {
GVertex *v0, *v1;
std::list<GFace *> l_faces;
+ std::set<GFace *> bl_faces;
+
public:
GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1);
virtual ~GEdge();
@@ -58,10 +62,16 @@ class GEdge : public GEntity {
// get the point for the given parameter location
virtual GPoint point(double p) const = 0;
-
+
// true if the edge contains the given parameter
virtual bool containsParam(double pt) const;
+ // get the position for the given parameter location
+ virtual SVector3 position(double p) const {
+ GPoint gp = point(p);
+ return SVector3(gp.x(),gp.y(),gp.z());
+ }
+
// get first derivative of edge at the given parameter
virtual SVector3 firstDer(double par) const = 0;
@@ -119,6 +129,24 @@ class GEdge : public GEntity {
// true if entity is periodic in the "dim" direction.
virtual bool periodic(int dim) const { return v0 == v1; }
+ // true if edge is used in hyperbolic layer on face gf
+ virtual bool inHyperbolicLayer(GFace* gf) {
+ return bl_faces.find(gf) != bl_faces.end();
+ }
+
+ virtual void flagInHyperbolicLayer(GFace* gf) {bl_faces.insert(gf);}
+
+ // get bounds of parametric coordinate
+ virtual Range<double> parBounds(int i) const = 0;
+
+ // return the point on the face closest to the given point
+ virtual GPoint closestPoint(const SPoint3 & queryPoint,double& param) const;
+
+ // try to reparametrise the point on the edge
+ virtual double parFromPoint(const SVector3& Q) const;
+
+ virtual bool XYZToU(const SVector3& Q,double& t,double relax=0.5) const;
+
struct {
char Method;
double coeffTransfinite;
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index f5ca9880de..09e443a9bd 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -99,6 +99,9 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
if(part) m->getMeshPartitions().insert(part);
}
+
+// extern void printDistortion(GModel *m, const char *nm,double& minR,double& maxR);
+
int GModel::readMSH(const std::string &name)
{
FILE *fp = fopen(name.c_str(), "rb");
@@ -215,6 +218,7 @@ int GModel::readMSH(const std::string &name)
if(fread(uv, sizeof(double), 1, fp) != 1) return 0;
if(swap) SwapBytes((char*)uv, sizeof(double), 1);
}
+ Msg::Warning("Creating edge vertex \n");
newVertex = new MEdgeVertex(xyz[0], xyz[1], xyz[2], ge, uv[0], -1.0, num);
}
else if (iClasDim == 2){
@@ -386,6 +390,14 @@ int GModel::readMSH(const std::string &name)
storePhysicalTagsInEntities(this, i, physicals[i]);
fclose(fp);
+
+ char nm[256] = "distorzione.pos";
+
+ double minR;
+ double maxR;
+
+ // printDistortion(this,nm,minR,maxR);
+
return postpro ? 2 : 1;
}
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 46f399ed82..41770ebb2e 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -361,9 +361,17 @@ bool reparamMeshVertexOnEdge(const MVertex *v, const GEdge *ge, double ¶m)
param = bounds.low();
else if(ge->getEndVertex() && ge->getEndVertex()->mesh_vertices[0] == v)
param = bounds.high();
- else
- v->getParameter(0, param);
-
+ else
+ param = ge->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+
+ if (v->onWhat() == ge) {
+ const MEdgeVertex* ve = dynamic_cast<const MEdgeVertex*>(v);
+ if (!ve) param = ge->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+ else v->getParameter(0, param);
+ }
+ else
+ param = ge->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+
if(param < 1.e6) return true;
return false;
}
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 3873f65405..d2fcfb5596 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -10,6 +10,7 @@
#include <set>
#include "SPoint2.h"
#include "SPoint3.h"
+#include "GmshMessage.h"
class GEntity;
class GEdge;
@@ -80,7 +81,7 @@ class MVertex{
// get/set the parent entity
inline GEntity* onWhat() const { return _ge; }
- inline void setEntity(GEntity *ge) { _ge = ge; }
+ inline void setEntity(GEntity *ge) {_ge = ge; }
// get/set the number
inline int getNum() const { return _num; }
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 54dcf1fbfa..cb966f103d 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -32,3 +32,7 @@ SVector3 discreteEdge::firstDer(double par) const
return SVector3();
}
+Range<double> discreteEdge::parBounds(int i) const {
+ Msg::Error("Cannot specify bounds for parametric coordinate on discrete edge");
+ return Range<double>(0,0);
+}
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index 44ce0bf4cc..fc24c51275 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -16,6 +16,7 @@ class discreteEdge : public GEdge {
virtual GeomType geomType() const { return DiscreteCurve; }
virtual GPoint point(double p) const;
virtual SVector3 firstDer(double par) const;
+ virtual Range<double> parBounds(int) const;
};
#endif
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index be7f91e87e..907122278a 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -237,11 +237,15 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
else
ve.insert(ve.end(), edgeVertices[p].rbegin(), edgeVertices[p].rend());
}
- else{
- MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
+
+ else{
+
+ MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
+ std::vector<MVertex*> temp;
+
double u0 = 0., u1 = 0., US[100];
bool reparamOK = true;
- if(!linear){
+ if(!linear) {
reparamOK &= reparamMeshVertexOnEdge(v0, ge, u0);
if(ge->periodic(0) && v1 == ge->getEndVertex()->mesh_vertices[0])
u1 = ge->parBounds(0).high();
@@ -251,40 +255,43 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
double relax = 1.;
while (1){
if(computeEquidistantParameters(ge, u0, u1, nPts + 2, US, relax))
- break;
+ break;
relax /= 2.0;
if(relax < 1.e-2)
break;
}
if(relax < 1.e-2)
- Msg::Warning("Failed to compute equidistant parameters (relax = %g)",
- relax);
+ Msg::Warning("Failed to compute equidistant parameters (relax = %g) for edge %d-%d",
+ relax,v0->getNum(),v1->getNum());
}
}
- std::vector<MVertex*> temp;
for(int j = 0; j < nPts; j++){
const double t = (double)(j + 1)/(nPts + 1);
+
double uc = (1. - t) * u0 + t * u1; // can be wrong, that's ok
MVertex *v;
if(linear || !reparamOK || uc < u0 || uc > u1){
+ Msg::Warning("We don't have a valid parameter on curve %d-%d",v0->getNum(),v1->getNum());
// we don't have a (valid) parameter on the curve
SPoint3 pc = edge.interpolate(t);
v = new MVertex(pc.x(), pc.y(), pc.z(), ge);
}
- else {
+ else {
GPoint pc = ge->point(US[j + 1]);
- v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[j + 1]);
- if (displ2D){
- SPoint3 pc2 = edge.interpolate(t);
- displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
- displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
- }
+ v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[j + 1]);
+
+ if (displ2D){
+ SPoint3 pc2 = edge.interpolate(t);
+ displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ }
}
temp.push_back(v);
// this destroys the ordering of the mesh vertices on the edge
ge->mesh_vertices.push_back(v);
ve.push_back(v);
}
+
if(edge.getVertex(0) == edge.getMinVertex())
edgeVertices[p].insert(edgeVertices[p].end(), temp.begin(), temp.end());
else
@@ -1019,12 +1026,18 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// then create new second order vertices/elements
edgeContainer edgeVertices;
faceContainer faceVertices;
+
+ Msg::StatusBar(1, true, "Meshing edges order %d...", order);
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
setHighOrder(*it, edgeVertices, linear, nPts, displ2D, displ3D);
+
+ Msg::StatusBar(1, true, "Meshing faces %d...", order);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,
displ2D, displ3D);
+ Msg::StatusBar(1, true, "Finished meshing order %d...", order);
+
// now we smooth mesh the internal vertices of the faces
// we do that model face by model face
std::vector<MElement*> bad;
--
GitLab