From 7b72b2a6f867011bc8e4aff10ce65f0e4ebc6e52 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Wed, 19 Nov 2008 20:34:02 +0000
Subject: [PATCH] cleanup high order
---
Common/Makefile | 6 +-
Fltk/Makefile | 25 +-
Geo/MVertex.cpp | 136 ++++
Geo/MVertex.h | 6 +
Geo/Makefile | 36 +-
Graphics/Makefile | 17 +-
Mesh/HighOrder.cpp | 1233 +++++-----------------------------
Mesh/HighOrder.h | 15 +-
Mesh/Makefile | 217 +++---
Mesh/Refine.cpp | 33 +-
Mesh/gmshSmoothHighOrder.cpp | 631 ++++++++++++++++-
Mesh/gmshSmoothHighOrder.h | 2 +-
Mesh/meshGFace.cpp | 3 +-
Numeric/Makefile | 23 +-
Plugin/Makefile | 4 +-
Post/Makefile | 23 +-
16 files changed, 1149 insertions(+), 1261 deletions(-)
diff --git a/Common/Makefile b/Common/Makefile
index 7fe6b6d542..9409f140c1 100644
--- a/Common/Makefile
+++ b/Common/Makefile
@@ -163,9 +163,9 @@ Visibility.o: Visibility.cpp Visibility.h GmshDefines.h ../Geo/GVertex.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEdge.h ../Geo/GFace.h ../Geo/GRegion.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h ../Parser/Parser.h
Trackball.o: Trackball.cpp Trackball.h
VertexArray.o: VertexArray.cpp VertexArray.h ../Geo/SVector3.h \
diff --git a/Fltk/Makefile b/Fltk/Makefile
index 002ad0a074..fb7456437b 100644
--- a/Fltk/Makefile
+++ b/Fltk/Makefile
@@ -113,9 +113,9 @@ GUI_Projection.o: GUI_Projection.cpp ../Geo/GModelIO_Fourier.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h ../Graphics/Draw.h ../Common/Options.h \
../Post/ColorTable.h ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/PartitionOptions.h ../Common/StringUtils.h \
@@ -136,9 +136,9 @@ GUI_Classifier.o: GUI_Classifier.cpp GUI_Classifier.h ../Common/GmshUI.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h GUI.h Opengl_Window.h Colorbar_Window.h \
../Post/ColorTable.h Popup_Button.h SpherePosition_Widget.h \
../Mesh/Field.h ../Post/PView.h Shortcut_Window.h ../Numeric/Numeric.h \
@@ -160,9 +160,9 @@ Callbacks.o: Callbacks.cpp ../Common/GmshUI.h ../Common/GmshMessage.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
../Geo/GeoStringInterface.h ../Geo/findLinks.h ../Mesh/Generator.h \
../Mesh/HighOrder.h ../Graphics/Draw.h ../Graphics/SelectBuffer.h \
../Post/PView.h ../Post/PViewOptions.h ../Post/ColorTable.h \
@@ -201,9 +201,10 @@ Opengl_Window.o: Opengl_Window.cpp ../Common/GmshUI.h \
../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h GUI.h \
Opengl_Window.h Colorbar_Window.h ../Post/ColorTable.h Popup_Button.h \
SpherePosition_Widget.h ../Mesh/Field.h ../Post/PView.h \
- ../Geo/MElement.h ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h
+ ../Geo/MElement.h ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h \
+ ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshMatrix.h
Colorbar_Window.o: Colorbar_Window.cpp ../Common/GmshUI.h GUI.h \
Opengl_Window.h Colorbar_Window.h ../Post/ColorTable.h Popup_Button.h \
SpherePosition_Widget.h ../Mesh/Field.h ../Post/PView.h \
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index 7bda5e677d..33d0d15950 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -6,6 +6,7 @@
#include <string.h>
#include <math.h>
#include "MVertex.h"
+#include "GVertex.h"
#include "GEdge.h"
#include "GFace.h"
#include "GFaceCompound.h"
@@ -228,3 +229,138 @@ void parametricCoordinates(const MVertex *ver, const GFace *gf, double &u, doubl
v = p.y();
}
}
+
+static void getAllParameters(MVertex *v, GFace *gf, std::vector<SPoint2> ¶ms)
+{
+ params.clear();
+ if(v->onWhat()->dim() == 0){
+ GVertex *gv = (GVertex*)v->onWhat();
+ std::list<GEdge*> ed = gv->edges();
+ bool seam = false;
+ for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++){
+ if((*it)->isSeam(gf)) {
+ Range<double> range = (*it)->parBounds(0);
+ if (gv == (*it)->getBeginVertex()){
+ params.push_back((*it)->reparamOnFace(gf, range.low(),-1));
+ params.push_back((*it)->reparamOnFace(gf, range.low(), 1));
+ }
+ else if (gv == (*it)->getEndVertex()){
+ params.push_back((*it)->reparamOnFace(gf, range.high(),-1));
+ params.push_back((*it)->reparamOnFace(gf, range.high(), 1));
+ }
+ else{
+ Msg::Warning("Strange!");
+ }
+ seam = true;
+ }
+ }
+ if (!seam)
+ params.push_back(gv->reparamOnFace(gf, 1));
+ }
+ else if(v->onWhat()->dim() == 1){
+ GEdge *ge = (GEdge*)v->onWhat();
+ double UU;
+ v->getParameter(0, UU);
+ params.push_back(ge->reparamOnFace(gf, UU, 1));
+ if(ge->isSeam(gf))
+ params.push_back(ge->reparamOnFace(gf, UU, -1));
+ }
+ else{
+ double UU, VV;
+ if(v->onWhat() == gf && v->getParameter(0, UU) && v->getParameter(1, VV))
+ params.push_back(SPoint2(UU, VV));
+ }
+}
+
+bool reparamMeshVerticesOnFace(MVertex *v1, MVertex *v2, GFace *gf,
+ SPoint2 ¶m1, SPoint2 ¶m2)
+{
+ std::vector<SPoint2> p1, p2;
+ getAllParameters(v1, gf, p1);
+ getAllParameters(v2, gf, p2);
+ if (p1.size() == 1 && p2.size() == 1){
+ param1 = p1[0];
+ param2 = p2[0];
+ return true;
+ }
+ else if (p1.size() == 1 && p2.size() == 2){
+ double d1 = (p1[0].x() - p2[0].x())*(p1[0].x() - p2[0].x())+
+ (p1[0].x() - p2[0].y())*(p1[0].y() - p2[0].y());
+ double d2 = (p1[0].x() - p2[1].x())*(p1[0].x() - p2[1].x())+
+ (p1[0].x() - p2[1].y())*(p1[0].y() - p2[1].y());
+ param1 = p1[0];
+ param2 = d2 < d1 ? p2[1] : p2[0];
+ return true;
+ }
+ else if (p2.size() == 1 && p1.size() == 2){
+ double d1 = (p2[0].x() - p1[0].x())*(p2[0].x() - p1[0].x())+
+ (p2[0].x() - p1[0].y())*(p2[0].y() - p1[0].y());
+ double d2 = (p2[0].x() - p1[1].x())*(p2[0].x() - p1[1].x())+
+ (p2[0].x() - p1[1].y())*(p2[0].y() - p1[1].y());
+ param1 = d2 < d1 ? p1[1] : p1[0];
+ param2 = p2[0];
+ return true;
+ }
+ return false;
+}
+
+bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
+{
+ if (gf->geomType() == GEntity::CompoundSurface){
+ GFaceCompound *gfc = (GFaceCompound*) gf;
+ param = gfc->getCoordinates(v);
+ return true;
+ }
+
+ if(v->onWhat()->geomType() == GEntity::DiscreteCurve ||
+ v->onWhat()->geomType() == GEntity::BoundaryLayerCurve){
+ param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+ return true;
+ }
+
+ if(v->onWhat()->dim() == 0){
+ GVertex *gv = (GVertex*)v->onWhat();
+
+ // abort if we could be on a seam
+ std::list<GEdge*> ed = gv->edges();
+ for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++)
+ if((*it)->isSeam(gf)) return false;
+
+ param = gv->reparamOnFace(gf, 1);
+ }
+ else if(v->onWhat()->dim() == 1){
+ GEdge *ge = (GEdge*)v->onWhat();
+
+ // abort if we are on a seam (todo: try dir=-1 and compare)
+ if(ge->isSeam(gf))
+ return false;
+
+ double UU;
+ v->getParameter(0, UU);
+ param = ge->reparamOnFace(gf, UU, 1);
+ }
+ else{
+ double UU, VV;
+ if(v->onWhat() == gf && v->getParameter(0, UU) && v->getParameter(1, VV))
+ param = SPoint2(UU, VV);
+ else
+ return false;
+ // param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
+ }
+ return true;
+}
+
+bool reparamMeshVertexOnEdge(MVertex *v, GEdge *ge, double ¶m)
+{
+ param = 1.e6;
+ Range<double> bounds = ge->parBounds(0);
+ if(ge->getBeginVertex() && ge->getBeginVertex()->mesh_vertices[0] == v)
+ param = bounds.low();
+ else if(ge->getEndVertex() && ge->getEndVertex()->mesh_vertices[0] == v)
+ param = bounds.high();
+ else
+ v->getParameter(0, param);
+
+ if(param < 1.e6) return true;
+ return false;
+}
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index d058c7dab2..12d6373f00 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -8,9 +8,11 @@
#include <stdio.h>
#include <set>
+#include "SPoint2.h"
#include "SPoint3.h"
class GEntity;
+class GEdge;
class GFace;
class MVertex;
@@ -164,5 +166,9 @@ class MFaceVertex : public MVertex{
};
void parametricCoordinates(const MVertex *ver, const GFace *gf, double &u, double &v);
+bool reparamMeshVerticesOnFace(MVertex *v1, MVertex *v2, GFace *gf,
+ SPoint2 ¶m1, SPoint2 ¶m2);
+bool reparamMeshVertexOnFace(MVertex *v, GFace *gf, SPoint2 ¶m);
+bool reparamMeshVertexOnEdge(MVertex *v, GEdge *ge, double ¶m);
#endif
diff --git a/Geo/Makefile b/Geo/Makefile
index 50c0366964..137ca3ecbd 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -105,8 +105,8 @@ GFaceCompound.o: GFaceCompound.cpp GFaceCompound.h GFace.h GEntity.h \
../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h ../Geo/GRegion.h \
../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
- ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
../Numeric/FunctionSpace.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h ../Common/Octree.h \
../Common/OctreeInternals.h ../Numeric/gmshLinearSystemGmm.h \
@@ -285,41 +285,41 @@ findLinks.o: findLinks.cpp ../Common/GmshMessage.h GModel.h GVertex.h \
GEntity.h Range.h SPoint3.h SBoundingBox3d.h GPoint.h SPoint2.h GEdge.h \
SVector3.h GFace.h GEdgeLoop.h Pair.h GRegion.h ../Common/TreeUtils.h \
../Common/avl.h ../Common/ListUtils.h
-MVertex.o: MVertex.cpp MVertex.h SPoint3.h GEdge.h GEntity.h Range.h \
- SBoundingBox3d.h GVertex.h GPoint.h SPoint2.h SVector3.h GFace.h \
+MVertex.o: MVertex.cpp MVertex.h SPoint2.h SPoint3.h GVertex.h GEntity.h \
+ Range.h SBoundingBox3d.h GPoint.h GEdge.h SVector3.h GFace.h \
GEdgeLoop.h Pair.h GFaceCompound.h ../Common/GmshMessage.h \
../Common/StringUtils.h
GaussQuadratureTri.o: GaussQuadratureTri.cpp MElement.h \
- ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshDefines.h MVertex.h SPoint2.h SPoint3.h MEdge.h \
+ SVector3.h MFace.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h GaussLegendreSimplex.h
GaussQuadratureQuad.o: GaussQuadratureQuad.cpp MElement.h \
- ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshDefines.h MVertex.h SPoint2.h SPoint3.h MEdge.h \
+ SVector3.h MFace.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h GaussLegendreSimplex.h \
../Numeric/GaussLegendre1D.h
GaussQuadratureTet.o: GaussQuadratureTet.cpp MElement.h \
- ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshDefines.h MVertex.h SPoint2.h SPoint3.h MEdge.h \
+ SVector3.h MFace.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h GaussLegendreSimplex.h
GaussQuadratureHex.o: GaussQuadratureHex.cpp MElement.h \
- ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshDefines.h MVertex.h SPoint2.h SPoint3.h MEdge.h \
+ SVector3.h MFace.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h GaussLegendreSimplex.h \
../Numeric/GaussLegendre1D.h
GaussLegendreSimplex.o: GaussLegendreSimplex.cpp MElement.h \
- ../Common/GmshDefines.h MVertex.h SPoint3.h MEdge.h SVector3.h MFace.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshDefines.h MVertex.h SPoint2.h SPoint3.h MEdge.h \
+ SVector3.h MFace.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h GaussLegendreSimplex.h \
../Numeric/GaussLegendre1D.h
-MFace.o: MFace.cpp MFace.h MVertex.h SPoint3.h SVector3.h \
+MFace.o: MFace.cpp MFace.h MVertex.h SPoint2.h SPoint3.h SVector3.h \
../Numeric/Numeric.h ../Numeric/NumericEmbedded.h ../Common/Context.h \
../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
MElement.o: MElement.cpp MElement.h ../Common/GmshDefines.h MVertex.h \
- SPoint3.h MEdge.h SVector3.h MFace.h ../Common/GmshMessage.h \
+ SPoint2.h SPoint3.h MEdge.h SVector3.h MFace.h ../Common/GmshMessage.h \
../Numeric/FunctionSpace.h ../Common/GmshMatrix.h GEntity.h Range.h \
- SBoundingBox3d.h GFace.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h \
- SPoint2.h Pair.h ../Common/StringUtils.h ../Numeric/Numeric.h \
+ SBoundingBox3d.h GFace.h GPoint.h GEdgeLoop.h GEdge.h GVertex.h Pair.h \
+ ../Common/StringUtils.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h ../Numeric/GaussLegendre1D.h \
../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h \
../Mesh/qualityMeasures.h ../Mesh/meshGFaceDelaunayInsertion.h \
diff --git a/Graphics/Makefile b/Graphics/Makefile
index 7114f5e54b..c41eb1c988 100644
--- a/Graphics/Makefile
+++ b/Graphics/Makefile
@@ -79,12 +79,13 @@ Mesh.o: Mesh.cpp ../Common/GmshMessage.h ../Common/GmshUI.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h Draw.h \
- ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h \
- ../Common/OS.h gl2ps.h ../Common/VertexArray.h ../Common/SmoothData.h \
- ../Post/PView.h ../Post/PViewData.h
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
+ Draw.h ../Common/Context.h ../Geo/CGNSOptions.h \
+ ../Mesh/PartitionOptions.h ../Common/OS.h gl2ps.h \
+ ../Common/VertexArray.h ../Common/SmoothData.h ../Post/PView.h \
+ ../Post/PViewData.h
Geom.o: Geom.cpp ../Common/GmshUI.h Draw.h ../Geo/SBoundingBox3d.h \
../Geo/SPoint3.h ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/PartitionOptions.h gl2ps.h ../Common/VertexArray.h \
@@ -113,8 +114,8 @@ SelectBuffer.o: SelectBuffer.cpp ../Common/GmshMessage.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Geo/MVertex.h \
- ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
../Numeric/FunctionSpace.h ../Common/GmshMatrix.h Draw.h \
../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h \
SelectBuffer.h ../Common/VertexArray.h
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index be5de5123d..c3042e4b0c 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -7,201 +7,56 @@
// Koen Hillewaert
//
+#include <iostream>
#include "HighOrder.h"
#include "gmshSmoothHighOrder.h"
-#include "gmshLaplace.h"
-#include "gmshLinearSystemGmm.h"
-#include "gmshAssembler.h"
-#include "meshGFaceOptimize.h"
#include "MElement.h"
#include "GmshMessage.h"
#include "OS.h"
#include "Numeric.h"
#include "Context.h"
-#include "GFaceCompound.h"
#include "GmshMatrix.h"
#include "FunctionSpace.h"
-#include <iostream>
#define SQU(a) ((a)*(a))
extern Context_T CTX;
-// for each pair of vertices (an edge), we build a list of vertices
-// that are the high order representation of the edge. The ordering of
-// vertices in the list is supposed to be (by construction) consistent
-// with the ordering of the pair.
-typedef std::map<std::pair<MVertex*,MVertex*>, std::vector<MVertex*> > edgeContainer;
-
-// for each face (a list of vertices) we build a list of vertices that
-// are the high order representation of the face
-// typedef std::map<std::vector<MVertex*>, std::vector<MVertex*>> faceContainer;
-typedef std::map<MFace, std::vector<MVertex*>,Less_Face> faceContainer;
-
-bool mappingIsInvertible (MTetrahedron *e)
+static bool mappingIsInvertible(MTetrahedron *e)
{
- if (e->getPolynomialOrder() == 1)return 1.0;
+ if (e->getPolynomialOrder() == 1) return 1.0;
double mat[3][3];
- e->getPrimaryJacobian(0.,0.,0., mat);
+ e->getPrimaryJacobian(0., 0., 0., mat);
double det0 = det3x3(mat);
IntPt *pts;
int npts;
-
- e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
+ e->getIntegrationPoints(e->getPolynomialOrder(), &npts, &pts);
- for (int i=0;i<npts;i++){
-
+ for (int i = 0; i < npts; i++){
const double u = pts[i].pt[0];
const double v = pts[i].pt[1];
const double w = pts[i].pt[2];
-
e->getJacobian(u, v, w, mat);
-
double detN = det3x3(mat);
-
if (det0 * detN <= 0.) return false;
}
- const Double_Matrix& points = e->getFunctionSpace()->points;
-
- for (int i=0;i<e->getNumPrimaryVertices();i++) {
+ const Double_Matrix &points = e->getFunctionSpace()->points;
+ for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
const double u = points(i,0);
const double v = points(i,1);
const double w = points(i,2);
-
- e->getJacobian(u,v,w,mat);
-
+ e->getJacobian(u, v, w, mat);
double detN = det3x3(mat);
-
if (det0 * detN <= 0.) return false;
}
return true;
}
-static void getAllParameters (MVertex *v, GFace *gf, std::vector<SPoint2> ¶ms){
- params.clear();
- if(v->onWhat()->dim() == 0){
- GVertex *gv = (GVertex*)v->onWhat();
- std::list<GEdge*> ed = gv->edges();
- bool seam = false;
- for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++)
- if((*it)->isSeam(gf)) {
- Range<double> range = (*it)->parBounds(0);
- if (gv == (*it)->getBeginVertex()){
- params.push_back((*it)->reparamOnFace(gf, range.low(),-1));
- params.push_back((*it)->reparamOnFace(gf, range.low(), 1));
- }
- else if (gv == (*it)->getEndVertex()){
- params.push_back((*it)->reparamOnFace(gf, range.high(),-1));
- params.push_back((*it)->reparamOnFace(gf, range.high(), 1));
- }
- else{
- printf("Strange !\n");
- }
- seam = true;
- }
- if (!seam)
- params.push_back(gv->reparamOnFace(gf, 1));
- }
- else if(v->onWhat()->dim() == 1){
- GEdge *ge = (GEdge*)v->onWhat();
- double UU;
- v->getParameter(0, UU);
- params.push_back(ge->reparamOnFace(gf, UU, 1));
- if(ge->isSeam(gf))
- params.push_back(ge->reparamOnFace(gf, UU, -1));
- }
- else{
- double UU, VV;
- if(v->onWhat() == gf && v->getParameter(0, UU) && v->getParameter(1, VV))
- params.push_back(SPoint2(UU, VV));
- }
-}
-
-bool reparamOnFace(MVertex *v1, MVertex *v2, GFace *gf, SPoint2 ¶m1, SPoint2 ¶m2){
- std::vector<SPoint2> p1,p2;
- getAllParameters (v1,gf,p1);
- getAllParameters (v2,gf,p2);
- if (p1.size() == 1 && p2.size() == 1){
- param1 = p1[0];
- param2 = p2[0];
- return true;
- }
- else if (p1.size() == 1 && p2.size() == 2){
- double d1 = (p1[0].x() - p2[0].x())*(p1[0].x() - p2[0].x())+
- (p1[0].x() - p2[0].y())*(p1[0].y() - p2[0].y());
- double d2 = (p1[0].x() - p2[1].x())*(p1[0].x() - p2[1].x())+
- (p1[0].x() - p2[1].y())*(p1[0].y() - p2[1].y());
- param1 = p1[0];
- param2 = d2 < d1 ? p2[1] : p2[0];
- return true;
- }
- else if (p2.size() == 1 && p1.size() == 2){
- double d1 = (p2[0].x() - p1[0].x())*(p2[0].x() - p1[0].x())+
- (p2[0].x() - p1[0].y())*(p2[0].y() - p1[0].y());
- double d2 = (p2[0].x() - p1[1].x())*(p2[0].x() - p1[1].x())+
- (p2[0].x() - p1[1].y())*(p2[0].y() - p1[1].y());
- param1 = d2 < d1 ? p1[1] : p1[0];
- param2 = p2[0];
- return true;
- }
- //printf("%d %d \n",p1.size(),p2.size());
- return false;
-}
-
-// FIXME: this should NOT be in HighOrder.cpp...
-bool reparamOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
-{
-
- if (gf->geomType() == GEntity::CompoundSurface){
- GFaceCompound *gfc = (GFaceCompound*) gf;
- param = gfc->getCoordinates(v);
- return true;
- }
-
- if(v->onWhat()->geomType() == GEntity::DiscreteCurve ||
- v->onWhat()->geomType() == GEntity::BoundaryLayerCurve){
- param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
- return true;
- }
-
- if(v->onWhat()->dim() == 0){
- GVertex *gv = (GVertex*)v->onWhat();
-
- // abort if we could be on a seam
-
- std::list<GEdge*> ed = gv->edges();
- for(std::list<GEdge*>::iterator it = ed.begin(); it != ed.end(); it++)
- if((*it)->isSeam(gf)) return false;
-
- param = gv->reparamOnFace(gf, 1);
- }
- else if(v->onWhat()->dim() == 1){
- GEdge *ge = (GEdge*)v->onWhat();
-
- // abort if we are on a seam (todo: try dir=-1 and compare)
- if(ge->isSeam(gf))
- return false;
-
- double UU;
- v->getParameter(0, UU);
- param = ge->reparamOnFace(gf, UU, 1);
- }
- else{
- double UU, VV;
- if(v->onWhat() == gf && v->getParameter(0, UU) && v->getParameter(1, VV))
- param = SPoint2(UU, VV);
- else
- return false;
- // param = gf->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
- }
- return true;
-}
-
// The aim here is to build a polynomial representation that consist
// in polynomial segments of equal length
@@ -217,7 +72,8 @@ static void myresid(int N, GEdge *ge, double *u, Double_Vector &r)
for (int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
}
-bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N, double *u, double underRelax)
+static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
+ double *u, double underRelax)
{
const double PRECISION = 1.e-6;
const int MAX_ITER = 50;
@@ -264,21 +120,15 @@ bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N, double
for (int i = 0; i < M; i++){
u[i+1] -= underRelax*DU(i);
}
- // printf("N %d M %d u1 = %g u0 = %g uN1 = %22.15E uN = %22.15E\n",
- // N, M, u[1], u0, u[N - 1], uN);
if (u[1] < u0) break;
if (u[N - 2] > uN) break;
double newt_norm = DU.norm();
- // printf("%22.15E\n",newt_norm);
- if (newt_norm < PRECISION) { /*printf("ok %g\n",underRelax);*/ return true; }
+ if (newt_norm < PRECISION) {
+ return true;
+ }
}
- // FAILED, use equidistant in param space
- // printf("failed %g\n",underRelax);
- // for (int i = 1; i < N; i++){
- // u[i] = u[i - 1] + du;
- // }
return false;
}
@@ -294,15 +144,15 @@ static void myresid(int N, GFace *gf, double *u, double *v, Double_Vector &r)
for (int i = 0; i < N - 2; i++) r(i) = L[i + 1] - L[i];
}
-bool computeEquidistantParameters(GFace *gf, double u0, double uN, double v0, double vN,
- int N, double *u, double *v)
+static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
+ double v0, double vN, int N,
+ double *u, double *v)
{
const double PRECISION = 1.e-6;
const int MAX_ITER = 50;
const double eps = 1.e-4;
double t[100];
-
// initialize the points by equal subdivision of geodesics
u[0] = u0;
v[0] = v0;
@@ -363,34 +213,19 @@ bool computeEquidistantParameters(GFace *gf, double u0, double uN, double v0, do
if (newt_norm < PRECISION) return true;
}
// FAILED, use equidistant in param space
- for (int i = 1; i < N; i++){
- t[i] = (double)i / (N - 1);
- SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i]);
- u[i] = p.x();
- v[i] = p.y();
- }
- return false;
-}
-
-bool reparamOnEdge(MVertex *v, GEdge *ge, double ¶m)
-{
- param = 1.e6;
- Range<double> bounds = ge->parBounds(0);
- if(ge->getBeginVertex() && ge->getBeginVertex()->mesh_vertices[0] == v)
- param = bounds.low();
- else if(ge->getEndVertex() && ge->getEndVertex()->mesh_vertices[0] == v)
- param = bounds.high();
- else
- v->getParameter(0, param);
-
- if(param < 1.e6) return true;
+ for (int i = 1; i < N; i++){
+ t[i] = (double)i / (N - 1);
+ SPoint2 p = gf->geodesic(SPoint2(u0, v0), SPoint2(uN, vN), t[i]);
+ u[i] = p.x();
+ v[i] = p.y();
+ }
return false;
}
-void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
- edgeContainer &edgeVertices, bool linear, int nPts = 1,
- gmshHighOrderSmoother *displ2D = 0,
- gmshHighOrderSmoother *displ3D = 0)
+static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
+ edgeContainer &edgeVertices, bool linear,
+ int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0)
{
for(int i = 0; i < ele->getNumEdges(); i++){
MEdge edge = ele->getEdge(i);
@@ -403,29 +238,31 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
}
else{
MVertex *v0 = edge.getVertex(0), *v1 = edge.getVertex(1);
-
double u0 = 0., u1 = 0.;
bool reparamOK = true;
if(!linear && ge->geomType() != GEntity::DiscreteCurve &&
ge->geomType() != GEntity::BoundaryLayerCurve){
- reparamOK &= reparamOnEdge(v0, ge, u0);
-
- if (ge->periodic(0) && v1 == ge->getEndVertex()->mesh_vertices[0]){
+ reparamOK &= reparamMeshVertexOnEdge(v0, ge, u0);
+ if (ge->periodic(0) && v1 == ge->getEndVertex()->mesh_vertices[0]){
Range<double> par = ge->parBounds(0);
u1 = par.high();
}
else
- reparamOK &= reparamOnEdge(v1, ge, u1);
+ reparamOK &= reparamMeshVertexOnEdge(v1, ge, u1);
}
double US[100];
if(reparamOK && !linear && ge->geomType() != GEntity::DiscreteCurve){
double relax = 1.;
while (1){
- if (computeEquidistantParameters(ge, u0, u1, nPts + 2, US,relax))break;
+ if(computeEquidistantParameters(ge, u0, u1, nPts + 2, US, relax))
+ break;
relax /= 2.0;
- if (relax < 1.e-2)break;
+ if (relax < 1.e-2)
+ break;
}
- if (relax < 1.e-2)Msg::Warning("failure in computing equidistant parameters %g",relax);
+ if (relax < 1.e-2)
+ Msg::Warning("Failed to compute equidistant parameters (relax = %g)",
+ relax);
}
std::vector<MVertex*> temp;
for(int j = 0; j < nPts; j++){
@@ -436,15 +273,15 @@ 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);
+ v->setParameter(0, t);
}
else {
- GPoint pc = ge->point(US[j+1]);
- v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[j+1]);
+ 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()));
+ displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
}
}
temp.push_back(v);
@@ -459,10 +296,10 @@ void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
}
}
-void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
- edgeContainer &edgeVertices, bool linear, int nPts = 1,
- gmshHighOrderSmoother *displ2D = 0,
- gmshHighOrderSmoother *displ3D = 0)
+static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
+ edgeContainer &edgeVertices, bool linear,
+ int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0)
{
for(int i = 0; i < ele->getNumEdges(); i++){
MEdge edge = ele->getEdge(i);
@@ -480,13 +317,10 @@ void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
if(!linear &&
gf->geomType() != GEntity::DiscreteSurface &&
gf->geomType() != GEntity::BoundaryLayerSurface){
-
- reparamOK = reparamOnFace(v0, v1, gf, p0, p1);
- // reparamOK &= reparamOnFace(v0, gf, p0);
- // reparamOK &= reparamOnFace(v1, gf, p1);
+ reparamOK = reparamMeshVerticesOnFace(v0, v1, gf, p0, p1);
}
double US[100], VS[100];
- if(reparamOK && !linear && gf->geomType() != GEntity::DiscreteCurve){
+ if(reparamOK && !linear && gf->geomType() != GEntity::DiscreteSurface){
computeEquidistantParameters(gf, p0[0], p1[0], p0[1], p1[1], nPts + 2, US, VS);
}
std::vector<MVertex*> temp;
@@ -517,14 +351,10 @@ void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
}
}
-void getEdgeVertices(GRegion *gr, MElement *ele,
- std::vector<MVertex*> &ve,
- std::set<MVertex*>& blocked,
- edgeContainer &edgeVertices,
- bool linear,
- int nPts = 1,
- gmshHighOrderSmoother *displ2D = 0,
- gmshHighOrderSmoother *displ3D = 0)
+static void getEdgeVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &ve,
+ std::set<MVertex*> &blocked, edgeContainer &edgeVertices,
+ bool linear, int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0)
{
for(int i = 0; i < ele->getNumEdges(); i++){
MEdge edge = ele->getEdge(i);
@@ -534,7 +364,7 @@ void getEdgeVertices(GRegion *gr, MElement *ele,
ve.insert(ve.end(), edgeVertices[p].begin(), edgeVertices[p].end());
else
ve.insert(ve.end(), edgeVertices[p].rbegin(), edgeVertices[p].rend());
- blocked.insert(edgeVertices[p].begin(),edgeVertices[p].end());
+ blocked.insert(edgeVertices[p].begin(), edgeVertices[p].end());
blocked.insert(edge.getMinVertex());
blocked.insert(edge.getMaxVertex());
}
@@ -556,18 +386,13 @@ void getEdgeVertices(GRegion *gr, MElement *ele,
}
}
-void getFaceVertices(GFace *gf,
- MElement *incomplete,
- MElement *ele,
- std::vector<MVertex*> &vf,
- faceContainer &faceVertices,
- bool linear, int nPts = 1,
- gmshHighOrderSmoother *displ2D = 0,
- gmshHighOrderSmoother *displ3D = 0) {
+static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
+ std::vector<MVertex*> &vf, faceContainer &faceVertices,
+ bool linear, int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0) {
Double_Matrix points;
int start = 0;
-
switch (nPts){
case 2 :
@@ -583,21 +408,14 @@ void getFaceVertices(GFace *gf,
start = 15;
break;
default :
- // do nothing (e.g. for quad faces)
+ // do nothing (e.g. for 2nd order tri faces or for quad faces)
break;
}
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
- //std::vector<MVertex*> p;
- // face.getOrderedVertices(p);
- //if(faceVertices.count(p)){
- // vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
- // }
-
faceContainer::iterator fIter = faceVertices.find(face);
-
- if(fIter!=faceVertices.end()){
+ if(fIter != faceVertices.end()){
vf.insert(vf.end(), fIter->second.begin(), fIter->second.end());
}
else{
@@ -606,14 +424,12 @@ void getFaceVertices(GFace *gf,
if(!linear &&
gf->geomType() != GEntity::DiscreteSurface &&
gf->geomType() != GEntity::BoundaryLayerSurface){
- for (int k=0;k<incomplete->getNumVertices(); k++){
- reparamOK &= reparamOnFace(incomplete->getVertex(k), gf, pts[k]);
+ for (int k = 0; k < incomplete->getNumVertices(); k++){
+ reparamOK &= reparamMeshVertexOnFace(incomplete->getVertex(k), gf, pts[k]);
}
}
if(face.getNumVertices() == 3){ // triangles
-
- std::vector<MVertex*>& vtcs = faceVertices[face];
-
+ std::vector<MVertex*> &vtcs = faceVertices[face];
for(int k = start ; k < points.size1() ; k++){
MVertex *v;
const double t1 = points(k, 0);
@@ -623,22 +439,12 @@ void getFaceVertices(GFace *gf,
v = new MVertex(pc.x(), pc.y(), pc.z(), gf);
}
else{
-
- // SPoint3 pos;
- // incomplete->pnt(t1,t2,0.,pos);
- // double X(pos.x());
- // double Y(pos.y());
- // double Z(pos.z());
-
-
- double X(0),Y(0),Z(0),GUESS[2]={0,0};
- double sf[256] ; incomplete->getShapeFunctions(t1,t2,0,sf);
- for (int j=0; j<incomplete->getNumVertices(); j++){
-
-
+ double X(0), Y(0), Z(0), GUESS[2] = {0, 0};
+ double sf[256] ;
+ incomplete->getShapeFunctions(t1, t2, 0, sf);
+ for (int j = 0; j < incomplete->getNumVertices(); j++){
MVertex *vt = incomplete->getVertex(j);
-
- X += sf[j] * vt->x();
+ X += sf[j] * vt->x();
Y += sf[j] * vt->y();
Z += sf[j] * vt->z();
if (reparamOK){
@@ -646,34 +452,31 @@ void getFaceVertices(GFace *gf,
GUESS[1] += sf[j] * pts[j][1];
}
}
- if (reparamOK){
- GPoint gp = gf->closestPoint(SPoint3(X,Y,Z),GUESS);
+ if(reparamOK){
+ GPoint gp = gf->closestPoint(SPoint3(X, Y, Z), GUESS);
if (gp.g()){
v = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
}
else{
- v = new MVertex(X,Y,Z, gf);
+ v = new MVertex(X, Y, Z, gf);
}
}
else{
- v = new MVertex(X,Y,Z, gf);
+ v = new MVertex(X, Y, Z, gf);
}
if (displ3D){
SPoint3 pc2 = face.interpolate(t1, t2);
- displ3D->add(v,SVector3(pc2.x(),pc2.y(),pc2.z()));
+ displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
}
}
// should be expensive -> induces a new search each time
- // faceVertices[p].push_back(v);
vtcs.push_back(v);
gf->mesh_vertices.push_back(v);
vf.push_back(v);
}
}
else if(face.getNumVertices() == 4){ // quadrangles
-
- std::vector<MVertex*>& vtcs = faceVertices[face];
-
+ std::vector<MVertex*> &vtcs = faceVertices[face];
for(int j = 0; j < nPts; j++){
for(int k = 0; k < nPts; k++){
MVertex *v;
@@ -696,7 +499,6 @@ void getFaceVertices(GFace *gf,
GPoint pc = gf->point(uc, vc);
v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
}
- // faceVertices[p].push_back(v);
vtcs.push_back(v);
gf->mesh_vertices.push_back(v);
vf.push_back(v);
@@ -707,11 +509,9 @@ void getFaceVertices(GFace *gf,
}
}
-
-void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
- faceContainer &faceVertices, bool linear, int nPts = 1)
+static void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
+ faceContainer &faceVertices, bool linear, int nPts = 1)
{
-
Double_Matrix points;
int start = 0;
@@ -729,37 +529,28 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
start = 15;
break;
default :
- // do nothing (e.g. for quad faces)
+ // do nothing (e.g. for 2nd order tri faces or for quad faces)
break;
}
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
- // std::vector<MVertex*> p;
-// face.getOrderedVertices(p);
-// if(faceVertices.count(p)){
-// vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
-// }
-
faceContainer::iterator fIter = faceVertices.find(face);
-
- if (fIter!=faceVertices.end()) {
+ if (fIter != faceVertices.end()) {
vf.insert(vf.end(), fIter->second.begin(), fIter->second.end());
}
else{
-
- std::vector<MVertex*>& vtcs = faceVertices[face];
-
+ std::vector<MVertex*> &vtcs = faceVertices[face];
SPoint2 p0, p1, p2, p3;
bool reparamOK = true;
if(!linear &&
gf->geomType() != GEntity::DiscreteSurface &&
gf->geomType() != GEntity::BoundaryLayerSurface){
- reparamOK &= reparamOnFace(ele->getVertex(0), gf, p0);
- reparamOK &= reparamOnFace(ele->getVertex(1), gf, p1);
- reparamOK &= reparamOnFace(ele->getVertex(2), gf, p2);
+ reparamOK &= reparamMeshVertexOnFace(ele->getVertex(0), gf, p0);
+ reparamOK &= reparamMeshVertexOnFace(ele->getVertex(1), gf, p1);
+ reparamOK &= reparamMeshVertexOnFace(ele->getVertex(2), gf, p2);
if(face.getNumVertices() == 4)
- reparamOK &= reparamOnFace(ele->getVertex(3), gf, p3);
+ reparamOK &= reparamMeshVertexOnFace(ele->getVertex(3), gf, p3);
}
if(face.getNumVertices() == 3){ // triangles
for(int k = start ; k < points.size1() ; k++){
@@ -778,7 +569,6 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
v->setParameter (0,uc);
v->setParameter (1,vc);
}
- // faceVertices[p].push_back(v);
vtcs.push_back(v);
gf->mesh_vertices.push_back(v);
vf.push_back(v);
@@ -807,7 +597,6 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
GPoint pc = gf->point(uc, vc);
v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
}
- // faceVertices[p].push_back(v);
vtcs.push_back(v);
gf->mesh_vertices.push_back(v);
vf.push_back(v);
@@ -818,57 +607,39 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
}
}
-void reorientTrianglePoints(std::vector<MVertex*>& vtcs,int orientation,bool swap) {
-
+static void reorientTrianglePoints(std::vector<MVertex*> &vtcs, int orientation,
+ bool swap)
+{
if (vtcs.size() == 1) return;
size_t nbPts = vtcs.size();
- if (nbPts > 3) Msg::Error("Interior face nodes reorientation not supported for order > 4");
+ if (nbPts > 3)
+ Msg::Error("Interior face nodes reorientation not supported for order > 4");
+
std::vector<MVertex*> tmp(nbPts);
// rotation
- // --------
-
// --- interior "principal vertices"
-
- // for (int i=0;i<3;i++) tmp[(i+orientation)%3] = vtcs[i];
-
- for (int i=0;i<3;i++) tmp[(i+orientation)%3] = vtcs[i];
-
-
+ for (int i = 0; i < 3; i++) tmp[(i + orientation) % 3] = vtcs[i];
+
// normal swap
- // -----------
-
if (swap) {
// --- interior "principal vertices"
-
-
vtcs[orientation] = tmp[orientation];
- vtcs[(orientation + 1) % 3] = tmp[(orientation+2)%3];
- vtcs[(orientation + 2) % 3] = tmp[(orientation+1)%3];
-
- // for (int i=0;i<2;i++) for (int j=0;j<2-i;j++) vtcs[i*2+j] = tmp[i+j*2];
+ vtcs[(orientation + 1) % 3] = tmp[(orientation + 2) % 3];
+ vtcs[(orientation + 2) % 3] = tmp[(orientation + 1) % 3];
}
-
// no swap
- // -------
-
else vtcs = tmp;
}
-
// KH: check face orientation wrt element ...
-void getFaceVertices(GRegion *gr, MElement *ele,
- std::vector<MVertex*> &vf,
- std::set<MVertex*>& blocked,
- faceContainer &faceVertices,
- edgeContainer& edgeVertices,
- bool linear, int nPts = 1)
+static void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
+ std::set<MVertex*> &blocked, faceContainer &faceVertices,
+ edgeContainer &edgeVertices, bool linear, int nPts = 1)
{
-
-
Double_Matrix points;
int start = 0;
@@ -887,70 +658,56 @@ void getFaceVertices(GRegion *gr, MElement *ele,
start = 15;
break;
default :
- // do nothing (e.g. for quad faces)
+ // do nothing (e.g. for 2nd order tri faces or for quad faces)
break;
}
-
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
-
faceContainer::iterator fIter = faceVertices.find(face);
-
- if (fIter!=faceVertices.end()) {
-
+ if (fIter != faceVertices.end()) {
int orientation;
bool swap;
-
std::vector<MVertex*> vtcs = fIter->second;
- if (fIter->first.computeCorrespondence(face,orientation,swap)) {
- reorientTrianglePoints(vtcs,orientation,swap);
+ if (fIter->first.computeCorrespondence(face, orientation, swap)) {
+ reorientTrianglePoints(vtcs, orientation, swap);
}
- else Msg::Error("Error in face lookup for recuperation of high order face nodes");
+ else
+ Msg::Error("Error in face lookup for recuperation of high order face nodes");
vf.insert(vf.end(), vtcs.begin(), vtcs.end());
-
blocked.insert(vtcs.begin(),vtcs.end());
blocked.insert(face.getVertex(0));
blocked.insert(face.getVertex(1));
blocked.insert(face.getVertex(2));
}
else{
- std::vector<MVertex*>& vtcs = faceVertices[face];
-
+ std::vector<MVertex*> &vtcs = faceVertices[face];
if(face.getNumVertices() == 3 && nPts > 1){ // triangles
-
- // construct incomplete element to take into account curved edges on surface boundaries
-
+ // construct incomplete element to take into account curved
+ // edges on surface boundaries
std::vector<MVertex*> hoEdgeNodes;
-
- for (int i=0;i<3;i++) {
-
+ for (int i = 0; i < 3; i++) {
MVertex* v0 = face.getVertex(i);
- MVertex* v1 = face.getVertex((i+1)%3);
-
- edgeContainer::iterator eIter = edgeVertices.find(std::pair<MVertex*,MVertex*>(std::min(v0,v1),std::max(v0,v1)));
-
- if (eIter == edgeVertices.end()) {
- printf("Could not find ho nodes for an edge");
- }
-
- if (v0 == eIter->first.first) hoEdgeNodes.insert(hoEdgeNodes.end(),eIter->second.begin(),eIter->second.end());
- else hoEdgeNodes.insert(hoEdgeNodes.end(),eIter->second.rbegin(),eIter->second.rend());
+ MVertex* v1 = face.getVertex((i + 1) % 3);
+ edgeContainer::iterator eIter = edgeVertices.find
+ (std::pair<MVertex*,MVertex*>(std::min(v0, v1), std::max(v0, v1)));
+ if (eIter == edgeVertices.end())
+ Msg::Error("Could not find ho nodes for an edge");
+ if (v0 == eIter->first.first)
+ hoEdgeNodes.insert(hoEdgeNodes.end(), eIter->second.begin(),
+ eIter->second.end());
+ else
+ hoEdgeNodes.insert(hoEdgeNodes.end(), eIter->second.rbegin(),
+ eIter->second.rend());
}
-
- MTriangleN incomplete(face.getVertex(0),
- face.getVertex(1),
- face.getVertex(2),
- hoEdgeNodes,nPts+1);
-
- for (int k=start;k<points.size1();k++) {
-
- double t1 = points(k,0);
- double t2 = points(k,1);
-
+ MTriangleN incomplete(face.getVertex(0), face.getVertex(1),
+ face.getVertex(2), hoEdgeNodes, nPts + 1);
+ for (int k = start; k < points.size1(); k++) {
+ double t1 = points(k, 0);
+ double t2 = points(k, 1);
SPoint3 pos;
- incomplete.pnt(t1,t2,0,pos);
- MVertex* v = new MVertex(pos.x(),pos.y(),pos.z(),gr);
+ incomplete.pnt(t1, t2, 0, pos);
+ MVertex* v = new MVertex(pos.x(), pos.y(), pos.z(), gr);
vtcs.push_back(v);
gr->mesh_vertices.push_back(v);
vf.push_back(v);
@@ -964,7 +721,6 @@ void getFaceVertices(GRegion *gr, MElement *ele,
double t2 = 2. * (double)(k + 1) / (nPts + 1) - 1.;
SPoint3 pc = face.interpolate(t1, t2);
MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
- // faceVertices[p].push_back(v);
vtcs.push_back(v);
gr->mesh_vertices.push_back(v);
vf.push_back(v);
@@ -975,14 +731,9 @@ void getFaceVertices(GRegion *gr, MElement *ele,
}
}
-
-void getRegionVertices(GRegion *gr,
- MElement *incomplete,
- MElement *ele,
- std::vector<MVertex*> &vr,
- bool linear, int nPts = 1)
+static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
+ std::vector<MVertex*> &vr, bool linear, int nPts = 1)
{
-
Double_Matrix points;
int start = 0;
@@ -996,8 +747,8 @@ void getRegionVertices(GRegion *gr,
start = 52;
break;
default :
+ // done: return!
return;
- break;
}
for(int k = start ; k < points.size1() ; k++){
@@ -1005,20 +756,17 @@ void getRegionVertices(GRegion *gr,
const double t1 = points(k, 0);
const double t2 = points(k, 1);
const double t3 = points(k, 2);
-
SPoint3 pos;
incomplete->pnt(t1,t2,t3,pos);
v = new MVertex(pos.x(),pos.y(),pos.z(),gr);
-
gr->mesh_vertices.push_back(v);
vr.push_back(v);
}
}
-
-void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
- int nbPts = 1, gmshHighOrderSmoother *displ2D = 0,
- gmshHighOrderSmoother *displ3D = 0)
+static void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
+ int nbPts = 1, gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0)
{
std::vector<MLine*> lines2;
for(unsigned int i = 0; i < ge->lines.size(); i++){
@@ -1035,10 +783,10 @@ void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
ge->deleteVertexArrays();
}
-void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
- faceContainer &faceVertices, bool linear, bool incomplete,
- int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
- gmshHighOrderSmoother *displ3D = 0)
+static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
+ faceContainer &faceVertices, bool linear, bool incomplete,
+ int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0)
{
std::vector<MTriangle*> triangles2;
for(unsigned int i = 0; i < gf->triangles.size(); i++){
@@ -1057,14 +805,9 @@ void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
ve, nPts + 1));
}
else{
- if (0 && gf->geomType() == GEntity::Plane){
- getFaceVertices(gf, t, vf, faceVertices, linear, nPts);
- }
- else{
- MTriangleN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2),
- ve, nPts + 1);
- getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts,displ2D,displ3D);
- }
+ MTriangleN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ ve, nPts + 1);
+ getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts,displ2D,displ3D);
ve.insert(ve.end(), vf.begin(), vf.end());
triangles2.push_back
(new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
@@ -1097,21 +840,19 @@ void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
gf->deleteVertexArrays();
}
-void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
- faceContainer &faceVertices, bool linear, bool incomplete,
- int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
- gmshHighOrderSmoother *displ3D = 0)
+static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
+ faceContainer &faceVertices, bool linear, bool incomplete,
+ int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0)
{
int nbCorr = 0;
std::vector<MTetrahedron*> tetrahedra2;
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
MTetrahedron *t = gr->tetrahedra[i];
-
std::set<MVertex*> blocked;
-
- std::vector<MVertex*> ve,vf,vr;
- getEdgeVertices(gr, t, ve, blocked,edgeVertices, linear, nPts,displ2D, displ3D);
+ std::vector<MVertex*> ve, vf, vr;
+ getEdgeVertices(gr, t, ve, blocked, edgeVertices, linear, nPts, displ2D, displ3D);
if (nPts == 1){
tetrahedra2.push_back
(new MTetrahedron10(t->getVertex(0), t->getVertex(1), t->getVertex(2),
@@ -1120,16 +861,14 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
else{
getFaceVertices(gr, t, vf, blocked, faceVertices, edgeVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
- MTetrahedronN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
- ve, nPts + 1);
+ MTetrahedronN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
+ ve, nPts + 1);
getRegionVertices(gr, &incpl, t, vr, linear, nPts);
ve.insert(ve.end(), vr.begin(), vr.end());
-
- MTetrahedron* n = new MTetrahedronN(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),ve, nPts + 1);
-
- if (!mappingIsInvertible(n)) {
+ MTetrahedron* n = new MTetrahedronN(t->getVertex(0), t->getVertex(1),
+ t->getVertex(2), t->getVertex(3), ve, nPts + 1);
+ if (!mappingIsInvertible(n))
Msg::Warning("Found invalid curved volume element (# %d in list) ",i);
- }
tetrahedra2.push_back (n);
}
delete t;
@@ -1137,21 +876,17 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
gr->tetrahedra = tetrahedra2;
std::vector<int> invalid;
-
if (nbCorr != 0) {
- for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
if (!mappingIsInvertible(gr->tetrahedra[i])) invalid.push_back(i);
- }
-
- if (invalid.size()!=0) {
- Msg::Warning("We have %d invalid elements remaining\n",(int) invalid.size());
+ if (invalid.size()) {
+ Msg::Warning("We have %d invalid elements remaining", (int)invalid.size());
std::vector<int>::iterator iIter = invalid.begin();
- for (;iIter!=invalid.end();++iIter) printf("%d ",*iIter);
- printf("\n");
+ for (; iIter != invalid.end(); ++iIter)
+ Msg::Warning("%d", *iIter);
}
}
-
std::vector<MHexahedron*> hexahedra2;
for(unsigned int i = 0; i < gr->hexahedra.size(); i++){
MHexahedron *h = gr->hexahedra[i];
@@ -1232,7 +967,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
}
template<class T>
-void setFirstOrder(GEntity *e, std::vector<T*> &elements)
+static void setFirstOrder(GEntity *e, std::vector<T*> &elements)
{
std::vector<T*> elements1;
for(unsigned int i = 0; i < elements.size(); i++){
@@ -1250,7 +985,7 @@ void setFirstOrder(GEntity *e, std::vector<T*> &elements)
e->deleteVertexArrays();
}
-void removeHighOrderVertices(GEntity *e)
+static void removeHighOrderVertices(GEntity *e)
{
std::vector<MVertex*> v1;
for(unsigned int i = 0; i < e->mesh_vertices.size(); i++){
@@ -1289,593 +1024,7 @@ void SetOrder1(GModel *m)
removeHighOrderVertices(*it);
}
-bool straightLine(std::vector<MVertex*> &l, MVertex *n1, MVertex *n2)
-{
- // x = a * t + b
- // x1 = b
- // x2 = a + b
- for(unsigned int i = 0; i < l.size(); i++){
- MVertex *v = l[i];
- double b = n1->x();
- double a = n2->x() - b;
- double t = (v->x() - b) / a;
- double by = n1->y();
- double ay = n2->y() - by;
- double y = ay * t + by;
- if(fabs(y-v->y()) > 1.e-07 * CTX.lc){
- return false;
- }
- }
- return true;
-}
-
-extern double mesh_functional_distorsion(MTriangle *t, double u, double v);
-
-void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
-{
- double mat[2][3];
- int n = 3;
- t->getPrimaryJacobian(0, 0, 0, mat);
- //t->getJacobian(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);
- double nn = sqrt(SQU(normal1[0]) + SQU(normal1[1]) + SQU(normal1[2]));
- for(int i = 0; i < n; i++){
- for(int k = 0; k < n - i; k++){
- t->getJacobian((double)i / (n - 1), (double)k / (n - 1), 0, mat);
- 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.) / nn;
- minJ = std::min(1. / det, std::min(det, minJ));
- maxJ = std::max(det, maxJ);
- }
- }
-}
-
-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 (unsigned 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 (unsigned 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(unsigned int i = 0; i < svd->ts.size(); i++){
- getMinMaxJac (svd->ts[i], minJ, maxJ);
- }
- for(unsigned 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 (unsigned int i = 0; i < svd->v.size(); i++){
- uv[i] += SMALL;
- dF[i] = (smooth_obj_HighOrderN(uv, data) - F) * LARGE;
- uv[i] -= SMALL;
- }
-}
-
-void optimizeNodeLocations(GFace *gf, smoothVertexDataHON &vdN, double eps = .2)
-{
- if(!vdN.v.size()) return;
- double uv[20];
- for (unsigned int i = 0; i < vdN.v.size(); i++){
- if (!vdN.v[i]->getParameter(0, uv[2 * i])){
- Msg::Error("Node location optimization failed");
- return;
- }
- if (!vdN.v[i]->getParameter(1, uv[2 * i + 1])){
- Msg::Error("Node location optimization failed");
- return;
- }
- }
-
- double F = -smooth_obj_HighOrderN(uv, &vdN);
- if (F < eps){
- double val;
- minimize_N(2 * vdN.v.size(),
- smooth_obj_HighOrderN,
- deriv_smoothing_objective_function_HighOrderN,
- &vdN, 1, uv,val);
- double Fafter = -smooth_obj_HighOrderN(uv, &vdN);
- printf("%12.5E %12.5E\n", F, Fafter);
- if (F < Fafter){
- for (unsigned 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();
- }
- }
- }
-}
-
-void optimizeHighOrderMeshInternalNodes(GFace *gf)
-{
- for(unsigned int i = 0; i < gf->triangles.size(); i++){
- MTriangle *t = gf->triangles[i];
- smoothVertexDataHON vdN;
- int start = t->getNumVertices() - t->getNumFaceVertices();
- for (int j=start;j<t->getNumVertices();j++)
- vdN.v.push_back(t->getVertex(j));
- vdN.gf = gf;
- vdN.ts.push_back(t);
- optimizeNodeLocations(gf, vdN, .9);
- }
-}
-
-bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices)
-{
- v2t_cont adjv;
- buildVertexToTriangle(gf->triangles, adjv);
-
- 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);
- }
- }
- /*
- 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);
-
- smoothVertexDataHON vdN;
- vdN.ts = it->second;
- for (int i=0;i<vdN.ts.size();i++){
- MTriangle *t = vdN.ts[i];
- }
-
- vdN.v = e;
- vdN.gf = gf;
-
- 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;
- }
- */
- 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() < 5){
- smoothVertexDataHON vdN;
- vdN.v = e;
- vdN.gf = gf;
- vdN.ts.clear();
- vdN.ts.push_back(t1);
- vdN.ts.push_back(t2);
- optimizeNodeLocations(gf, vdN);
- }
- }
- }
-
- return success;
-}
-
-double angle3Points ( MVertex *p1, MVertex *p2, MVertex *p3 ){
- SVector3 a(p1->x()-p2->x(),p1->y()-p2->y(),p1->z()-p2->z());
- SVector3 b(p3->x()-p2->x(),p3->y()-p2->y(),p3->z()-p2->z());
- SVector3 c = crossprod(a,b);
- double sinA = c.norm();
- double cosA = dot(a,b);
- // printf("%d %d %d -> %g %g\n",p1->iD,p2->iD,p3->iD,cosA,sinA);
- return atan2 (sinA,cosA);
-}
-
-/*
-A curvilinear edge smooth and swap
-
-*/
-
-typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MElement*> > edge2tris;
-
-void localHarmonicMapping(GModel *gm,
- MTriangle *t1 ,
- MTriangle *t2,
- MVertex *n1,
- MVertex *n2,
- MVertex *n3,
- MVertex *n4,
-// SPoint2 &np1,
-// SPoint2 &np2,
-// SPoint2 &np3,
-// SPoint2 &np4,
- std::vector<MVertex*> &e1,
- std::vector<MVertex*> &e2,
- std::vector<MVertex*> &e3,
- std::vector<MVertex*> &e4,
-// std::vector<SPoint2> &ep1,
-// std::vector<SPoint2> &ep2,
-// std::vector<SPoint2> &ep3,
-// std::vector<SPoint2> &ep4
- std::vector<MVertex*> &e) {
-
- gmshLinearSystemGmm *lsys = new gmshLinearSystemGmm;
- gmshAssembler myAssembler(lsys);
- gmshLaplaceTerm Laplace (gm,1.0,0);
-
- myAssembler.fixVertex ( n1 , 0 , 0 , -1.0);
- myAssembler.fixVertex ( n2 , 0 , 0 , -1.0);
- myAssembler.fixVertex ( n3 , 0 , 0 , 1.0);
- myAssembler.fixVertex ( n4 , 0 , 0 , 1.0);
- for (unsigned int i = 0; i < e1.size(); i++) myAssembler.fixVertex(e1[i], 0, 0, -1.0);
- for (unsigned int i = 0; i < e3.size(); i++) myAssembler.fixVertex(e3[i], 0, 0, 1.0);
- Laplace.addToMatrix(myAssembler,t1);
- Laplace.addToMatrix(myAssembler,t2);
- lsys->systemSolve();
-
- gmshLinearSystemGmm *lsys1 = new gmshLinearSystemGmm;
- gmshAssembler myAssembler1(lsys1);
- gmshLaplaceTerm Laplace1 (gm,1.0,1);
-
- myAssembler1.fixVertex ( n2 , 0 , 1 , -1.0);
- myAssembler1.fixVertex ( n3 , 0 , 1 , -1.0);
- myAssembler1.fixVertex ( n4 , 0 , 1 , 1.0);
- myAssembler1.fixVertex ( n1 , 0 , 1 , 1.0);
- for (unsigned int i = 0; i < e2.size(); i++) myAssembler1.fixVertex(e2[i], 0, 1, -1.0);
- for (unsigned int i = 0; i < e4.size(); i++) myAssembler1.fixVertex(e4[i], 0, 1, 1.0);
- Laplace1.addToMatrix(myAssembler1,t1);
- Laplace1.addToMatrix(myAssembler1,t2);
- lsys1->systemSolve();
-
- // now we have the stable high order harmonic mapping
- // we have to find points locations of vertices in e
- // that have coordinates (\xi, \xi)
-
- // this can be done by evaluating the
-
- for (unsigned int i = 0; i < e.size(); i++){
- MVertex *v = e[i];
- const double U = myAssembler.getDofValue (v, 0 ,0);
- const double V = myAssembler1.getDofValue (v, 0 ,1);
- printf("point %g %g -> %g %g\n",v->x(),v->y(),U,V);
- // we are in t1
- if (U >= V){
- const double ut = U;
- }
- }
-
-
- delete lsys ;
- delete lsys1;
-}
-
-
-
-void getParametricCoordnates ( GFace *gf,
- std::vector<MVertex*> &e,
- std::vector<SPoint2> ¶m){
- param.clear();
- for (unsigned int i = 0; i < e.size(); i++){
- double U,V;
- parametricCoordinates(e[i] , gf, U, V);
- param.push_back(SPoint2(U,V));
- }
-}
-
-static void curvilinearEdgeSwap (GFace *gf,
- // int nPts,
- edgeContainer &edgeVertices,
- edge2tris::iterator &it,
- edge2tris &e2t)
-{
- std::pair<MVertex*, MVertex*> edge = it->first;
- 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];
- MVertex *n1 = t1->getOtherVertex(n2, n4);
- MVertex *n3 = t2->getOtherVertex(n2, n4);
- std::vector<MVertex*> e1 = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n1, n2),std::max(n1, n2))];
- std::vector<MVertex*> e2 = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n2, n3),std::max(n2, n3))];
- std::vector<MVertex*> e3 = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n3, n4),std::max(n3, n4))];
- std::vector<MVertex*> e4 = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n4, n1),std::max(n4, n1))];
- std::vector<MVertex*> e = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n2, n4),std::max(n2, n4))];
- // std::vector<MVertex*> enew;
- // MLine temp (n1,n3);
- // should not add the nodes n the GFace here
- // getEdgeVertices(gf,&temp, enew, false, nPts);
- // get the parametric coordinates of the
- std::vector<SPoint2> ep1; getParametricCoordnates (gf,e1,ep1);
- std::vector<SPoint2> ep2; getParametricCoordnates (gf,e2,ep2);
- std::vector<SPoint2> ep3; getParametricCoordnates (gf,e3,ep3);
- std::vector<SPoint2> ep4; getParametricCoordnates (gf,e4,ep4);
- std::vector<SPoint2> ep; getParametricCoordnates (gf,e ,ep );
- // std::vector<SPoint2> epnew; getParametricCoordnates (gf,enew,epnew);
- localHarmonicMapping(gf->model(),t1,t2,n1,n2,n3,n4,e1,e2,e3,e4,e);
- }
-}
-
-bool smoothInternalEdgesb(GFace *gf, edgeContainer &edgeVertices)
-{
- 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);
- }
- }
-
- for(edge2tris::iterator it = e2t.begin(); it != e2t.end(); ++it){
- curvilinearEdgeSwap (gf,edgeVertices,it,e2t);
- }
- return true;
-}
-
-bool smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
-{
- 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;
-
- 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;
- 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];
- MVertex *n1 = t1->getOtherVertex(n2, n4);
- MVertex *n3 = t2->getOtherVertex(n2, n4);
-
- double ang1 = angle3Points(n1,n2,n3);
- double ang2 = angle3Points(n2,n3,n4);
- double ang3 = angle3Points(n3,n4,n1);
- double ang4 = angle3Points(n4,n1,n2);
- const double angleLimit = 3*M_PI/4.;
-
- if (ang1 < angleLimit && ang2 < angleLimit && ang3 < angleLimit && ang4 < angleLimit){
- if(n1 < n2)
- e1 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n1, n2)];
- else
- e1 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n1)];
- if(n2 < n3)
- e2 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n3)];
- else
- e2 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n3, n2)];
- if(n3 < n4)
- e3 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n3, n4)];
- else
- e3 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n3)];
- if(n4 < n1)
- e4 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n1)];
- else
- e4 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n1, n4)];
- if(n2 < n4)
- e = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n4)];
- else
- e = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n2)];
-
- 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++){
- 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 = %g min %g max %g\n",relax,minJloc,maxJloc);
-
- if (minJloc > minJ){
- kopt = k;
- minJ = minJloc;
- }
- }
- // kopt = 1;
- 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, std::vector<MElement*> & bad, double & minJGlob)
+static void checkHighOrderTriangles(GModel *m, std::vector<MElement*> &bad, double &minJGlob)
{
bad.clear();
minJGlob = 1.0;
@@ -1883,23 +1032,22 @@ void checkHighOrderTriangles(GModel *m, std::vector<MElement*> & bad, double & m
for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
MTriangle *t = (*it)->triangles[i];
double disto = t->distoShapeMeasure();
- minJGlob = std::min(minJGlob,disto);
+ minJGlob = std::min(minJGlob, disto);
if (disto < 0) bad.push_back(t);
}
}
if (minJGlob > 0) Msg::Info("Worst Element Smoothness %12.5E", minJGlob);
else Msg::Warning("Worst Element Smoothness %12.5E", minJGlob);
-}
+}
-void printJacobians(GModel *m, const char *nm)
+extern double mesh_functional_distorsion(MTriangle *t, double u, double v);
+
+static void printJacobians(GModel *m, const char *nm)
{
- return;
+ return;
const int n = 25;
- double D[n][n];
- double X[n][n];
- double Y[n][n];
- double Z[n][n];
+ double D[n][n], X[n][n], Y[n][n], Z[n][n];
FILE *f = fopen(nm,"w");
fprintf(f,"View \"\"{\n");
@@ -1916,40 +1064,25 @@ void printJacobians(GModel *m, const char *nm)
X[i][k] = pt.x();
Y[i][k] = pt.y();
Z[i][k] = pt.z();
- // printf("X[%d][%d] = %g %g\n",i,k,X[i][k],pt.x());
}
}
for(int i= 0; i < n -1; i++){
for(int k = 0; k < n - i -1; k++){
- fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
- X[i][k],Y[i][k],Z[i][k],
- X[i+1][k],Y[i+1][k],Z[i+1][k],
- X[i][k+1],Y[i][k+1],Z[i][k+1],
- D[i][k],
- D[i+1][k],
- D[i][k+1]);
-// fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
-// X[i][k],X[i+1][k],X[i][k+1],
-// Y[i][k],Y[i+1][k],Y[i][k+1],
-// Z[i][k],Z[i+1][k],Z[i][k+1],
-// D[i][k],
-// D[i+1][k],
-// D[i][k+1]);
+ fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
+ X[i][k],Y[i][k],Z[i][k],
+ X[i+1][k],Y[i+1][k],Z[i+1][k],
+ X[i][k+1],Y[i][k+1],Z[i][k+1],
+ D[i][k],
+ D[i+1][k],
+ D[i][k+1]);
if (i != n-2 && k != n - i -2)
- fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
- X[i+1][k],Y[i+1][k],Z[i+1][k],
- X[i+1][k+1],Y[i+1][k+1],Z[i+1][k+1],
- X[i][k+1],Y[i][k+1],Z[i][k+1],
- D[i+1][k],
- D[i+1][k+1],
- D[i][k+1]);
-// fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
-// X[i+1][k],X[i+1][k+1],X[i][k+1],
-// Y[i+1][k],Y[i+1][k+1],Y[i][k+1],
-// Z[i+1][k],Z[i+1][k+1],Z[i][k+1],
-// D[i+1][k],
-// D[i+1][k+1],
-// D[i][k+1]);
+ fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
+ X[i+1][k],Y[i+1][k],Z[i+1][k],
+ X[i+1][k+1],Y[i+1][k+1],Z[i+1][k+1],
+ X[i][k+1],Y[i][k+1],Z[i][k+1],
+ D[i+1][k],
+ D[i+1][k+1],
+ D[i][k+1]);
}
}
}
@@ -1982,7 +1115,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
int nPts = order - 1;
- Msg::StatusBar(1, true, "Generating High Order Nodes (q = %d) ...",order);
+ Msg::StatusBar(1, true, "Generating High Order Nodes (q = %d) ...", order);
double t1 = Cpu();
// first, make sure to remove any existsing second order vertices/elements
@@ -1999,9 +1132,10 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
edgeContainer edgeVertices;
faceContainer faceVertices;
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
- setHighOrder(*it, edgeVertices, linear, nPts,displ2D,displ3D);
+ setHighOrder(*it, edgeVertices, linear, nPts, displ2D, displ3D);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
- setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,displ2D, displ3D);
+ setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,
+ displ2D, displ3D);
// now we smooth mesh the internal vertices of the faces
@@ -2009,14 +1143,15 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
std::vector<MElement*> bad;
double worst;
if (displ2D){
- checkHighOrderTriangles(m,bad,worst);
+ checkHighOrderTriangles(m, bad, worst);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
- if ((*it)->geomType() == GEntity::Plane)displ2D->smooth(*it);
+ if ((*it)->geomType() == GEntity::Plane) displ2D->smooth(*it);
// will have to smooth in the planar coordinates, using the metric
}
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
- setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,displ2D, displ3D);
+ setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts,
+ displ2D, displ3D);
// smooth the 3D regions
if (displ3D){
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index 6e1d1a678d..54a9592f19 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -7,12 +7,21 @@
#define _HIGH_ORDER_H_
#include "GModel.h"
+#include "MFace.h"
+
+// for each pair of vertices (an edge), we build a list of vertices
+// that are the high order representation of the edge. The ordering of
+// vertices in the list is supposed to be (by construction) consistent
+// with the ordering of the pair.
+typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MVertex*> > edgeContainer;
+
+// for each face (a list of vertices) we build a list of vertices that
+// are the high order representation of the face
+// typedef std::map<std::vector<MVertex*>, std::vector<MVertex*>> faceContainer;
+typedef std::map<MFace, std::vector<MVertex*>, Less_Face> faceContainer;
void SetOrder1(GModel *m);
void SetOrderN(GModel *m, int order, bool linear=true, bool incomplete=false);
-void checkHighOrderTriangles(GModel *m);
-bool reparamOnFace(MVertex *v, GFace *gf, SPoint2 ¶m);
-
void Refine(GModel *m, bool linear=true);
#endif
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 62b220946e..d58491614b 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -87,10 +87,10 @@ Generator.o: Generator.cpp ../Common/GmshMessage.h ../Numeric/Numeric.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h meshGEdge.h \
- meshGFace.h meshGFaceBDS.h meshGRegion.h BackgroundMesh.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
+ meshGEdge.h meshGFace.h meshGFaceBDS.h meshGRegion.h BackgroundMesh.h \
BoundaryLayers.h HighOrder.h ../Post/PView.h ../Post/PViewData.h
Field.o: Field.cpp ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/PartitionOptions.h Field.h ../Post/PView.h ../Geo/SPoint3.h \
@@ -110,25 +110,29 @@ Field.o: Field.cpp ../Common/Context.h ../Geo/CGNSOptions.h \
../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
../Common/GmshMessage.h ../Post/OctreePost.h ../Common/Octree.h \
../Common/OctreeInternals.h ../Post/PViewDataList.h ../Post/PViewData.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h
-gmshSmoothHighOrder.o: gmshSmoothHighOrder.cpp gmshSmoothHighOrder.h \
- ../Geo/SVector3.h ../Geo/SPoint3.h ../Numeric/gmshAssembler.h \
- ../Numeric/gmshLinearSystem.h ../Numeric/gmshLaplace.h \
- ../Numeric/gmshTermOfFormulation.h ../Common/GmshMatrix.h \
- ../Common/Gmsh.h ../Common/GmshMessage.h ../Geo/GModel.h \
- ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
- ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
- ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/GFace.h \
- ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h \
- ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h \
- ../Geo/GEntity.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h
+gmshSmoothHighOrder.o: gmshSmoothHighOrder.cpp HighOrder.h \
+ ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+ ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+ ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
+ ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h \
+ ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h \
+ ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
+ ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+ ../Geo/SBoundingBox3d.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h meshGFaceOptimize.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
- ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h ../Numeric/gmshElasticity.h \
+ ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
+ ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshMatrix.h meshGFaceDelaunayInsertion.h \
+ gmshSmoothHighOrder.h ../Numeric/gmshAssembler.h \
+ ../Numeric/gmshLinearSystem.h ../Numeric/gmshLaplace.h \
+ ../Numeric/gmshTermOfFormulation.h ../Common/Gmsh.h \
+ ../Common/GmshMessage.h ../Numeric/gmshElasticity.h \
../Numeric/gmshTermOfFormulation.h ../Numeric/gmshLinearSystemGmm.h \
- ../Numeric/gmshLinearSystem.h
+ ../Numeric/gmshLinearSystem.h ../Common/Context.h ../Geo/CGNSOptions.h \
+ ../Mesh/PartitionOptions.h ../Numeric/Numeric.h \
+ ../Numeric/NumericEmbedded.h
meshGEdge.o: meshGEdge.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -138,12 +142,13 @@ meshGEdge.o: meshGEdge.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h meshGEdge.h ../Geo/MElement.h \
- ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
- ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Common/GmshMessage.h \
- ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h BackgroundMesh.h \
- ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h ../Common/Context.h \
- ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
+ ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+ ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshMatrix.h BackgroundMesh.h ../Numeric/Numeric.h \
+ ../Numeric/NumericEmbedded.h ../Common/Context.h ../Geo/CGNSOptions.h \
+ ../Mesh/PartitionOptions.h
meshGEdgeExtruded.o: meshGEdgeExtruded.cpp ../Geo/GModel.h \
../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -153,24 +158,25 @@ meshGEdgeExtruded.o: meshGEdgeExtruded.cpp ../Geo/GModel.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h
meshGFace.o: meshGFace.cpp meshGFace.h meshGFaceBDS.h \
meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
- ../Common/GmshMatrix.h meshGFaceQuadrilateralize.h meshGFaceOptimize.h \
- DivideAndConquer.h BackgroundMesh.h ../Geo/GVertex.h ../Geo/GEntity.h \
- ../Geo/Range.h ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h \
- ../Geo/SPoint3.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h \
- ../Geo/GEntity.h ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
- ../Geo/SPoint2.h ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h \
- ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
- ../Geo/Pair.h ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h \
- ../Geo/GFace.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Common/GmshMessage.h \
+ ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
+ meshGFaceQuadrilateralize.h meshGFaceOptimize.h DivideAndConquer.h \
+ BackgroundMesh.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+ ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
+ ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/GEdge.h ../Geo/GEntity.h \
+ ../Geo/GVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
+ ../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
+ ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
+ ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEdge.h ../Geo/GFace.h \
+ ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/PartitionOptions.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h BDS.h ../Post/PView.h qualityMeasures.h \
@@ -182,7 +188,7 @@ meshGFaceTransfinite.o: meshGFaceTransfinite.cpp meshGFace.h \
../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h \
../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MElement.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MElement.h \
../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/MEdge.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
@@ -198,16 +204,16 @@ meshGFaceExtruded.o: meshGFaceExtruded.cpp ../Geo/GModel.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h \
../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
meshGFaceBDS.o: meshGFaceBDS.cpp meshGFace.h meshGFaceOptimize.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
- ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h meshGFaceDelaunayInsertion.h BackgroundMesh.h \
../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -228,17 +234,17 @@ meshGFaceDelaunayInsertion.o: meshGFaceDelaunayInsertion.cpp BDS.h \
../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Common/GmshMessage.h \
BackgroundMesh.h meshGFaceDelaunayInsertion.h ../Geo/MElement.h \
- ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
- ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
- ../Common/GmshMatrix.h meshGFaceOptimize.h meshGFace.h \
- ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h ../Common/Context.h \
- ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
+ ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+ ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h meshGFaceOptimize.h \
+ meshGFace.h ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
+ ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h
meshGFaceOptimize.o: meshGFaceOptimize.cpp meshGFaceOptimize.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
- ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h meshGFaceDelaunayInsertion.h qualityMeasures.h \
../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -250,11 +256,11 @@ meshGFaceQuadrilateralize.o: meshGFaceQuadrilateralize.cpp \
meshGFaceQuadrilateralize.h ../Common/GmshMessage.h \
../Numeric/Numeric.h ../Numeric/NumericEmbedded.h \
meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
- meshGFaceOptimize.h meshGFaceBDS.h BDS.h ../Geo/GFace.h \
- ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshMatrix.h meshGFaceOptimize.h meshGFaceBDS.h BDS.h \
+ ../Geo/GFace.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/GEntity.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
@@ -262,9 +268,9 @@ meshGFaceQuadrilateralize.o: meshGFaceQuadrilateralize.cpp \
../Geo/Pair.h ../Post/PView.h
meshGRegion.o: meshGRegion.cpp meshGRegion.h \
meshGRegionDelaunayInsertion.h ../Geo/MElement.h \
- ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
- ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
- ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+ ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h BackgroundMesh.h qualityMeasures.h \
@@ -285,9 +291,9 @@ meshGRegion.o: meshGRegion.cpp meshGRegion.h \
meshGRegionDelaunayInsertion.o: meshGRegionDelaunayInsertion.cpp \
../Common/OS.h BackgroundMesh.h meshGRegion.h meshGRegionLocalMeshMod.h \
meshGRegionDelaunayInsertion.h ../Geo/MElement.h \
- ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
- ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/SPoint3.h \
- ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+ ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h qualityMeasures.h ../Geo/GModel.h \
@@ -305,7 +311,7 @@ meshGRegionTransfinite.o: meshGRegionTransfinite.cpp meshGFace.h \
../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MElement.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MElement.h \
../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/MEdge.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
@@ -320,9 +326,9 @@ meshGRegionExtruded.o: meshGRegionExtruded.cpp ../Geo/GModel.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h \
meshGFace.h meshGRegion.h ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/PartitionOptions.h
@@ -335,16 +341,16 @@ meshGRegionCarveHole.o: meshGRegionCarveHole.cpp ../Geo/GModel.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h
meshGRegionLocalMeshMod.o: meshGRegionLocalMeshMod.cpp \
meshGRegionLocalMeshMod.h meshGRegionDelaunayInsertion.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
- ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h BackgroundMesh.h qualityMeasures.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
@@ -372,7 +378,7 @@ qualityMeasures.o: qualityMeasures.cpp qualityMeasures.h BDS.h \
../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Common/GmshMessage.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MElement.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MElement.h \
../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/MEdge.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
@@ -386,9 +392,9 @@ BoundaryLayers.o: BoundaryLayers.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h BoundaryLayers.h ../Geo/ExtrudeParams.h \
../Common/SmoothData.h meshGEdge.h meshGFace.h
BDS.o: BDS.cpp ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h BDS.h \
@@ -399,9 +405,10 @@ BDS.o: BDS.cpp ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h BDS.h \
../Geo/SPoint3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
../Geo/SVector3.h ../Geo/Pair.h ../Post/PView.h ../Common/GmshMessage.h \
meshGFaceDelaunayInsertion.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h qualityMeasures.h
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
+ qualityMeasures.h
HighOrder.o: HighOrder.cpp HighOrder.h ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -410,18 +417,14 @@ HighOrder.o: HighOrder.cpp HighOrder.h ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h gmshSmoothHighOrder.h ../Numeric/gmshLaplace.h \
- ../Numeric/gmshTermOfFormulation.h ../Common/GmshMatrix.h \
- ../Common/Gmsh.h ../Common/GmshMessage.h ../Geo/MElement.h \
- ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
- ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
- ../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
- ../Numeric/gmshAssembler.h ../Numeric/gmshLinearSystem.h \
- meshGFaceOptimize.h meshGFaceDelaunayInsertion.h ../Common/OS.h \
- ../Numeric/Numeric.h ../Numeric/NumericEmbedded.h ../Common/Context.h \
- ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h ../Geo/GFaceCompound.h \
- ../Geo/GFace.h ../Geo/GEdge.h
+ ../Geo/SBoundingBox3d.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
+ gmshSmoothHighOrder.h ../Geo/MElement.h ../Common/GmshDefines.h \
+ ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/MFace.h ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshMatrix.h ../Common/OS.h ../Numeric/Numeric.h \
+ ../Numeric/NumericEmbedded.h ../Common/Context.h ../Geo/CGNSOptions.h \
+ ../Mesh/PartitionOptions.h
Partition.o: Partition.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -432,8 +435,8 @@ Partition.o: Partition.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h Partition.h PartitionObjects.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
- ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
../Common/GmshMatrix.h PartitionOptions.h
Refine.o: Refine.cpp HighOrder.h ../Geo/GModel.h ../Geo/GVertex.h \
@@ -444,11 +447,9 @@ Refine.o: Refine.cpp HighOrder.h ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h meshGFaceOptimize.h ../Geo/MElement.h \
- ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
- ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Common/GmshMessage.h \
- ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h \
- meshGFaceDelaunayInsertion.h ../Common/OS.h ../Numeric/Numeric.h \
- ../Numeric/NumericEmbedded.h ../Common/Context.h ../Geo/CGNSOptions.h \
- ../Mesh/PartitionOptions.h
+ ../Geo/SBoundingBox3d.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h ../Geo/MElement.h \
+ ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
+ ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshMatrix.h ../Common/OS.h
diff --git a/Mesh/Refine.cpp b/Mesh/Refine.cpp
index bd3cdae766..f453dfd57d 100644
--- a/Mesh/Refine.cpp
+++ b/Mesh/Refine.cpp
@@ -14,7 +14,7 @@
void Subdivide(GEdge *ge)
{
- std::vector<MLine *> lines2;
+ std::vector<MLine*> lines2;
for(unsigned int i = 0; i < ge->lines.size(); i++){
MLine *l = ge->lines[i];
if(l->getNumVertices() == 3){
@@ -164,13 +164,17 @@ void Subdivide(GRegion *gr)
if(p->getNumVertices() == 14){
// BASE
pyramids2.push_back
- (new MPyramid(p->getVertex(0), p->getVertex(5), p->getVertex(13), p->getVertex(6), p->getVertex(7)));
+ (new MPyramid(p->getVertex(0), p->getVertex(5), p->getVertex(13),
+ p->getVertex(6), p->getVertex(7)));
pyramids2.push_back
- (new MPyramid(p->getVertex(5), p->getVertex(1), p->getVertex(8), p->getVertex(13), p->getVertex(9)));
+ (new MPyramid(p->getVertex(5), p->getVertex(1), p->getVertex(8),
+ p->getVertex(13), p->getVertex(9)));
pyramids2.push_back
- (new MPyramid(p->getVertex(13), p->getVertex(8), p->getVertex(2), p->getVertex(10), p->getVertex(11)));
+ (new MPyramid(p->getVertex(13), p->getVertex(8), p->getVertex(2),
+ p->getVertex(10), p->getVertex(11)));
pyramids2.push_back
- (new MPyramid(p->getVertex(6), p->getVertex(13), p->getVertex(10), p->getVertex(3), p->getVertex(12)));
+ (new MPyramid(p->getVertex(6), p->getVertex(13), p->getVertex(10),
+ p->getVertex(3), p->getVertex(12)));
// Split remaining into tets
// Top
@@ -203,30 +207,19 @@ void Subdivide(GRegion *gr)
void Refine(GModel *m, bool linear)
{
- // Refine all edges in a mesh by inserting a point on the midpoint
- // then recreate edge, face, and region definitions
- //
- // If linear is set to true, new vertices are created by linear
- // interpolation between existing ones. If not, new vertices are
- // created on the exact geometry, provided that the geometrical
- // edges/faces are discretized with 1D/2D elements. (I.e., if there
- // are only 3D elements in the mesh then any new vertices on the
- // boundary will always be created by linear interpolation, whether
- // linear is set or not.)
-
Msg::StatusBar(1, true, "Generating refined mesh ...");
double t1 = Cpu();
- // first, set order to 2 to generate vertex positions
+ // Create 2nd order mesh (using "2nd order complete" elements) to
+ // generate vertex positions
SetOrderN(m, 2, linear, false);
- // Now subdivide entities to create linear mesh
+ // Subdivide the second order elements to create the refined linear
+ // mesh
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it)
Subdivide(*it);
-
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
Subdivide(*it);
-
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
Subdivide(*it);
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index c6d68fb56d..92d46d94c1 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -1,17 +1,33 @@
+// Gmsh - Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributor(s):
+// Koen Hillewaert
+//
+
+#include "HighOrder.h"
+#include "meshGFaceOptimize.h"
#include "gmshSmoothHighOrder.h"
-
#include "gmshAssembler.h"
#include "gmshLaplace.h"
#include "gmshElasticity.h"
#include "gmshLinearSystemGmm.h"
#include "GFace.h"
#include "GRegion.h"
+#include "Context.h"
+#include "Numeric.h"
+
+#define SQU(a) ((a)*(a))
+extern Context_T CTX;
-void getDistordedElements ( const std::vector<MElement*> & v,
- const double & threshold,
- std::vector<MElement*> & d,
- double &minD){
+void getDistordedElements(const std::vector<MElement*> & v,
+ const double & threshold,
+ std::vector<MElement*> & d,
+ double &minD)
+{
d.clear();
minD = 1;
for (unsigned int i = 0; i < v.size(); i++){
@@ -22,9 +38,10 @@ void getDistordedElements ( const std::vector<MElement*> & v,
}
}
-void addOneLayer ( const std::vector<MElement*> & v,
- std::vector<MElement*> & d ,
- std::vector<MElement*> & layer ){
+void addOneLayer(const std::vector<MElement*> & v,
+ std::vector<MElement*> & d ,
+ std::vector<MElement*> & layer )
+{
std::set<MVertex*> all;
for (unsigned int i = 0; i < d.size(); i++){
MElement *e = d[i];
@@ -55,8 +72,8 @@ void addOneLayer ( const std::vector<MElement*> & v,
}
}
-
-void gmshHighOrderSmoother::smooth ( GFace *gf) {
+void gmshHighOrderSmoother::smooth(GFace *gf)
+{
std::vector<MElement*> v;
v.insert(v.begin(), gf->triangles.begin(),gf->triangles.end());
@@ -65,7 +82,9 @@ void gmshHighOrderSmoother::smooth ( GFace *gf) {
v.size());
smooth(v);
}
-void gmshHighOrderSmoother::smooth ( GRegion *gr) {
+
+void gmshHighOrderSmoother::smooth(GRegion *gr)
+{
std::vector<MElement*> v;
v.insert(v.begin(), gr->tetrahedra.begin(),gr->tetrahedra.end());
v.insert(v.begin(), gr->hexahedra.begin(),gr->hexahedra.end());
@@ -75,9 +94,8 @@ void gmshHighOrderSmoother::smooth ( GRegion *gr) {
smooth(v);
}
-void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all) {
-
-
+void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
+{
gmshLinearSystemGmm *lsys = new gmshLinearSystemGmm;
gmshAssembler myAssembler(lsys);
gmshElasticityTerm El(0,1.0,.333,getTag());
@@ -196,5 +214,590 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all) {
}
+////////////////////////////////////////////////////////////////////////////////
+// OLD STUFF : STILL USED ?
+////////////////////////////////////////////////////////////////////////////////
+
+bool straightLine(std::vector<MVertex*> &l, MVertex *n1, MVertex *n2)
+{
+ // x = a * t + b
+ // x1 = b
+ // x2 = a + b
+ for(unsigned int i = 0; i < l.size(); i++){
+ MVertex *v = l[i];
+ double b = n1->x();
+ double a = n2->x() - b;
+ double t = (v->x() - b) / a;
+ double by = n1->y();
+ double ay = n2->y() - by;
+ double y = ay * t + by;
+ if(fabs(y-v->y()) > 1.e-07 * CTX.lc){
+ return false;
+ }
+ }
+ return true;
+}
+
+void getMinMaxJac (MTriangle *t, double &minJ, double &maxJ)
+{
+ double mat[2][3];
+ int n = 3;
+ t->getPrimaryJacobian(0, 0, 0, mat);
+ //t->getJacobian(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);
+ double nn = sqrt(SQU(normal1[0]) + SQU(normal1[1]) + SQU(normal1[2]));
+ for(int i = 0; i < n; i++){
+ for(int k = 0; k < n - i; k++){
+ t->getJacobian((double)i / (n - 1), (double)k / (n - 1), 0, mat);
+ 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.) / nn;
+ minJ = std::min(1. / det, std::min(det, minJ));
+ maxJ = std::max(det, maxJ);
+ }
+ }
+}
+
+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 (unsigned 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 (unsigned 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(unsigned int i = 0; i < svd->ts.size(); i++){
+ getMinMaxJac (svd->ts[i], minJ, maxJ);
+ }
+ for(unsigned 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 (unsigned int i = 0; i < svd->v.size(); i++){
+ uv[i] += SMALL;
+ dF[i] = (smooth_obj_HighOrderN(uv, data) - F) * LARGE;
+ uv[i] -= SMALL;
+ }
+}
+
+void optimizeNodeLocations(GFace *gf, smoothVertexDataHON &vdN, double eps = .2)
+{
+ if(!vdN.v.size()) return;
+ double uv[20];
+ for (unsigned int i = 0; i < vdN.v.size(); i++){
+ if (!vdN.v[i]->getParameter(0, uv[2 * i])){
+ Msg::Error("Node location optimization failed");
+ return;
+ }
+ if (!vdN.v[i]->getParameter(1, uv[2 * i + 1])){
+ Msg::Error("Node location optimization failed");
+ return;
+ }
+ }
+
+ double F = -smooth_obj_HighOrderN(uv, &vdN);
+ if (F < eps){
+ double val;
+ minimize_N(2 * vdN.v.size(),
+ smooth_obj_HighOrderN,
+ deriv_smoothing_objective_function_HighOrderN,
+ &vdN, 1, uv,val);
+ double Fafter = -smooth_obj_HighOrderN(uv, &vdN);
+ printf("%12.5E %12.5E\n", F, Fafter);
+ if (F < Fafter){
+ for (unsigned 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();
+ }
+ }
+ }
+}
+
+double angle3Points ( MVertex *p1, MVertex *p2, MVertex *p3 )
+{
+ SVector3 a(p1->x()-p2->x(),p1->y()-p2->y(),p1->z()-p2->z());
+ SVector3 b(p3->x()-p2->x(),p3->y()-p2->y(),p3->z()-p2->z());
+ SVector3 c = crossprod(a,b);
+ double sinA = c.norm();
+ double cosA = dot(a,b);
+ // printf("%d %d %d -> %g %g\n",p1->iD,p2->iD,p3->iD,cosA,sinA);
+ return atan2 (sinA,cosA);
+}
+
+/*
+A curvilinear edge smooth and swap
+*/
+
+typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MElement*> > edge2tris;
+
+void localHarmonicMapping(GModel *gm,
+ MTriangle *t1 ,
+ MTriangle *t2,
+ MVertex *n1,
+ MVertex *n2,
+ MVertex *n3,
+ MVertex *n4,
+// SPoint2 &np1,
+// SPoint2 &np2,
+// SPoint2 &np3,
+// SPoint2 &np4,
+ std::vector<MVertex*> &e1,
+ std::vector<MVertex*> &e2,
+ std::vector<MVertex*> &e3,
+ std::vector<MVertex*> &e4,
+// std::vector<SPoint2> &ep1,
+// std::vector<SPoint2> &ep2,
+// std::vector<SPoint2> &ep3,
+// std::vector<SPoint2> &ep4
+ std::vector<MVertex*> &e) {
+
+ gmshLinearSystemGmm *lsys = new gmshLinearSystemGmm;
+ gmshAssembler myAssembler(lsys);
+ gmshLaplaceTerm Laplace (gm,1.0,0);
+
+ myAssembler.fixVertex ( n1 , 0 , 0 , -1.0);
+ myAssembler.fixVertex ( n2 , 0 , 0 , -1.0);
+ myAssembler.fixVertex ( n3 , 0 , 0 , 1.0);
+ myAssembler.fixVertex ( n4 , 0 , 0 , 1.0);
+ for (unsigned int i = 0; i < e1.size(); i++) myAssembler.fixVertex(e1[i], 0, 0, -1.0);
+ for (unsigned int i = 0; i < e3.size(); i++) myAssembler.fixVertex(e3[i], 0, 0, 1.0);
+ Laplace.addToMatrix(myAssembler,t1);
+ Laplace.addToMatrix(myAssembler,t2);
+ lsys->systemSolve();
+
+ gmshLinearSystemGmm *lsys1 = new gmshLinearSystemGmm;
+ gmshAssembler myAssembler1(lsys1);
+ gmshLaplaceTerm Laplace1 (gm,1.0,1);
+
+ myAssembler1.fixVertex ( n2 , 0 , 1 , -1.0);
+ myAssembler1.fixVertex ( n3 , 0 , 1 , -1.0);
+ myAssembler1.fixVertex ( n4 , 0 , 1 , 1.0);
+ myAssembler1.fixVertex ( n1 , 0 , 1 , 1.0);
+ for (unsigned int i = 0; i < e2.size(); i++) myAssembler1.fixVertex(e2[i], 0, 1, -1.0);
+ for (unsigned int i = 0; i < e4.size(); i++) myAssembler1.fixVertex(e4[i], 0, 1, 1.0);
+ Laplace1.addToMatrix(myAssembler1,t1);
+ Laplace1.addToMatrix(myAssembler1,t2);
+ lsys1->systemSolve();
+
+ // now we have the stable high order harmonic mapping
+ // we have to find points locations of vertices in e
+ // that have coordinates (\xi, \xi)
+
+ // this can be done by evaluating the
+
+ for (unsigned int i = 0; i < e.size(); i++){
+ MVertex *v = e[i];
+ const double U = myAssembler.getDofValue (v, 0 ,0);
+ const double V = myAssembler1.getDofValue (v, 0 ,1);
+ printf("point %g %g -> %g %g\n",v->x(),v->y(),U,V);
+ // we are in t1
+ if (U >= V){
+ const double ut = U;
+ }
+ }
+
+
+ delete lsys ;
+ delete lsys1;
+}
+
+void optimizeHighOrderMeshInternalNodes(GFace *gf)
+{
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ MTriangle *t = gf->triangles[i];
+ smoothVertexDataHON vdN;
+ int start = t->getNumVertices() - t->getNumFaceVertices();
+ for (int j=start;j<t->getNumVertices();j++)
+ vdN.v.push_back(t->getVertex(j));
+ vdN.gf = gf;
+ vdN.ts.push_back(t);
+ optimizeNodeLocations(gf, vdN, .9);
+ }
+}
+
+bool optimizeHighOrderMesh(GFace *gf, edgeContainer &edgeVertices)
+{
+ v2t_cont adjv;
+ buildVertexToTriangle(gf->triangles, adjv);
+
+ 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);
+ }
+ }
+ /*
+ 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);
+
+ smoothVertexDataHON vdN;
+ vdN.ts = it->second;
+ for (int i=0;i<vdN.ts.size();i++){
+ MTriangle *t = vdN.ts[i];
+ }
+
+ vdN.v = e;
+ vdN.gf = gf;
+
+ 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;
+ }
+ */
+ 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() < 5){
+ smoothVertexDataHON vdN;
+ vdN.v = e;
+ vdN.gf = gf;
+ vdN.ts.clear();
+ vdN.ts.push_back(t1);
+ vdN.ts.push_back(t2);
+ optimizeNodeLocations(gf, vdN);
+ }
+ }
+ }
+
+ return success;
+}
+
+void getParametricCoordnates ( GFace *gf,
+ std::vector<MVertex*> &e,
+ std::vector<SPoint2> ¶m)
+{
+ param.clear();
+ for (unsigned int i = 0; i < e.size(); i++){
+ double U,V;
+ parametricCoordinates(e[i] , gf, U, V);
+ param.push_back(SPoint2(U,V));
+ }
+}
+
+static void curvilinearEdgeSwap (GFace *gf,
+ // int nPts,
+ edgeContainer &edgeVertices,
+ edge2tris::iterator &it,
+ edge2tris &e2t)
+{
+ std::pair<MVertex*, MVertex*> edge = it->first;
+ 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];
+ MVertex *n1 = t1->getOtherVertex(n2, n4);
+ MVertex *n3 = t2->getOtherVertex(n2, n4);
+ std::vector<MVertex*> e1 = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n1, n2),std::max(n1, n2))];
+ std::vector<MVertex*> e2 = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n2, n3),std::max(n2, n3))];
+ std::vector<MVertex*> e3 = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n3, n4),std::max(n3, n4))];
+ std::vector<MVertex*> e4 = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n4, n1),std::max(n4, n1))];
+ std::vector<MVertex*> e = edgeVertices[std::make_pair<MVertex*, MVertex*>(std::min(n2, n4),std::max(n2, n4))];
+ // std::vector<MVertex*> enew;
+ // MLine temp (n1,n3);
+ // should not add the nodes n the GFace here
+ // getEdgeVertices(gf,&temp, enew, false, nPts);
+ // get the parametric coordinates of the
+ std::vector<SPoint2> ep1; getParametricCoordnates (gf,e1,ep1);
+ std::vector<SPoint2> ep2; getParametricCoordnates (gf,e2,ep2);
+ std::vector<SPoint2> ep3; getParametricCoordnates (gf,e3,ep3);
+ std::vector<SPoint2> ep4; getParametricCoordnates (gf,e4,ep4);
+ std::vector<SPoint2> ep; getParametricCoordnates (gf,e ,ep );
+ // std::vector<SPoint2> epnew; getParametricCoordnates (gf,enew,epnew);
+ localHarmonicMapping(gf->model(),t1,t2,n1,n2,n3,n4,e1,e2,e3,e4,e);
+ }
+}
+
+bool smoothInternalEdgesb(GFace *gf, edgeContainer &edgeVertices)
+{
+ 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);
+ }
+ }
+
+ for(edge2tris::iterator it = e2t.begin(); it != e2t.end(); ++it){
+ curvilinearEdgeSwap (gf,edgeVertices,it,e2t);
+ }
+ return true;
+}
+
+bool smoothInternalEdges(GFace *gf, edgeContainer &edgeVertices)
+{
+ 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;
+
+ 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;
+ 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];
+ MVertex *n1 = t1->getOtherVertex(n2, n4);
+ MVertex *n3 = t2->getOtherVertex(n2, n4);
+
+ double ang1 = angle3Points(n1,n2,n3);
+ double ang2 = angle3Points(n2,n3,n4);
+ double ang3 = angle3Points(n3,n4,n1);
+ double ang4 = angle3Points(n4,n1,n2);
+ const double angleLimit = 3*M_PI/4.;
+
+ if (ang1 < angleLimit && ang2 < angleLimit && ang3 < angleLimit && ang4 < angleLimit){
+ if(n1 < n2)
+ e1 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n1, n2)];
+ else
+ e1 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n1)];
+ if(n2 < n3)
+ e2 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n3)];
+ else
+ e2 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n3, n2)];
+ if(n3 < n4)
+ e3 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n3, n4)];
+ else
+ e3 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n3)];
+ if(n4 < n1)
+ e4 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n1)];
+ else
+ e4 = edgeVertices[std::make_pair<MVertex*, MVertex*>(n1, n4)];
+ if(n2 < n4)
+ e = edgeVertices[std::make_pair<MVertex*, MVertex*>(n2, n4)];
+ else
+ e = edgeVertices[std::make_pair<MVertex*, MVertex*>(n4, n2)];
+
+ 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++){
+ 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 = %g min %g max %g\n",relax,minJloc,maxJloc);
+
+ if (minJloc > minJ){
+ kopt = k;
+ minJ = minJloc;
+ }
+ }
+ // kopt = 1;
+ 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;
+}
diff --git a/Mesh/gmshSmoothHighOrder.h b/Mesh/gmshSmoothHighOrder.h
index 5647cb0eb4..502852c264 100644
--- a/Mesh/gmshSmoothHighOrder.h
+++ b/Mesh/gmshSmoothHighOrder.h
@@ -9,6 +9,7 @@
#include <map>
#include <vector>
#include "SVector3.h"
+
class MVertex;
class MElement;
class GFace;
@@ -29,5 +30,4 @@ public:
int getTag() const {return _tag;}
};
-
#endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 24dd6c2f34..73248d70d3 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -389,7 +389,8 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool debug = true)
SBoundingBox3d bbox;
while(itv != all_vertices.end()){
MVertex *here = *itv;
- SPoint2 param ; reparamOnFace (here,gf,param);;
+ SPoint2 param;
+ reparamMeshVertexOnFace(here, gf, param);
U_[count] = param.x();
V_[count] = param.y();
(*itv)->setIndex(count);
diff --git a/Numeric/Makefile b/Numeric/Makefile
index 6f11d25daf..d7feb8c59c 100644
--- a/Numeric/Makefile
+++ b/Numeric/Makefile
@@ -54,8 +54,8 @@ Numeric.o: Numeric.cpp ../Common/GmshMessage.h Numeric.h \
NumericEmbedded.h
NumericEmbedded.o: NumericEmbedded.cpp NumericEmbedded.h \
../Common/GmshMessage.h
-gmshAssembler.o: gmshAssembler.cpp ../Geo/MVertex.h ../Geo/SPoint3.h \
- gmshAssembler.h gmshLinearSystem.h
+gmshAssembler.o: gmshAssembler.cpp ../Geo/MVertex.h ../Geo/SPoint2.h \
+ ../Geo/SPoint3.h gmshAssembler.h gmshLinearSystem.h
gmshTermOfFormulation.o: gmshTermOfFormulation.cpp \
gmshTermOfFormulation.h ../Common/GmshMatrix.h gmshFunction.h \
../Common/Gmsh.h ../Common/GmshMessage.h ../Geo/GModel.h \
@@ -67,9 +67,10 @@ gmshTermOfFormulation.o: gmshTermOfFormulation.cpp \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h gmshLinearSystem.h gmshAssembler.h
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h gmshLinearSystem.h \
+ gmshAssembler.h
gmshLaplace.o: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
../Common/GmshMatrix.h ../Common/Gmsh.h ../Common/GmshMessage.h \
../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
@@ -80,9 +81,9 @@ gmshLaplace.o: gmshLaplace.cpp gmshLaplace.h gmshTermOfFormulation.h \
../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h \
../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h
gmshElasticity.o: gmshElasticity.cpp gmshElasticity.h \
gmshTermOfFormulation.h ../Common/GmshMatrix.h ../Common/Gmsh.h \
../Common/GmshMessage.h ../Geo/GModel.h ../Geo/GVertex.h \
@@ -94,9 +95,9 @@ gmshElasticity.o: gmshElasticity.cpp gmshElasticity.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
- ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Numeric/FunctionSpace.h
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h \
+ ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Numeric/FunctionSpace.h
EigSolve.o: EigSolve.cpp
FunctionSpace.o: FunctionSpace.cpp FunctionSpace.h ../Common/GmshMatrix.h \
../Common/GmshDefines.h ../Common/GmshMessage.h
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 99e676226e..b7721bda26 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -170,8 +170,8 @@ Triangulate.o: Triangulate.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/discreteFace.h ../Geo/GModel.h \
../Geo/GFace.h ../Mesh/DivideAndConquer.h ../Common/GmshMessage.h \
- ../Geo/MVertex.h ../Geo/SPoint3.h Triangulate.h Plugin.h \
- ../Common/Options.h ../Post/ColorTable.h ../Post/PView.h \
+ ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h Triangulate.h \
+ Plugin.h ../Common/Options.h ../Post/ColorTable.h ../Post/PView.h \
../Post/PViewDataList.h ../Post/PViewData.h ../Common/ListUtils.h \
../Common/GmshMatrix.h ../Common/Context.h ../Geo/CGNSOptions.h \
../Mesh/PartitionOptions.h
diff --git a/Post/Makefile b/Post/Makefile
index 3bcfb016c1..d9f9eab52f 100644
--- a/Post/Makefile
+++ b/Post/Makefile
@@ -98,10 +98,11 @@ PViewDataGModel.o: PViewDataGModel.cpp PViewDataGModel.h PViewData.h \
../Geo/GPoint.h ../Geo/GEdgeLoop.h ../Geo/GEdge.h ../Geo/SPoint2.h \
../Geo/SVector3.h ../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h \
../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/MElement.h \
- ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint3.h \
- ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Common/GmshMessage.h \
- ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h ../Numeric/Numeric.h \
+ ../Common/GmshDefines.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+ ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Common/GmshMessage.h ../Numeric/FunctionSpace.h \
+ ../Common/GmshMatrix.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h
PViewDataGModelIO.o: PViewDataGModelIO.cpp ../Common/GmshMessage.h \
PViewDataGModel.h PViewData.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
@@ -112,11 +113,11 @@ PViewDataGModelIO.o: PViewDataGModelIO.cpp ../Common/GmshMessage.h \
../Geo/GFace.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/GEdgeLoop.h \
../Geo/GEdge.h ../Geo/SPoint2.h ../Geo/SVector3.h ../Geo/Pair.h \
../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.h ../Geo/MVertex.h ../Geo/SPoint3.h \
- ../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
- ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h \
- ../Geo/MVertex.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
- ../Common/GmshMatrix.h ../Numeric/Numeric.h \
+ ../Geo/SBoundingBox3d.h ../Geo/MVertex.h ../Geo/SPoint2.h \
+ ../Geo/SPoint3.h ../Geo/MElement.h ../Common/GmshDefines.h \
+ ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Numeric/FunctionSpace.h ../Common/GmshMatrix.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h ../Common/StringUtils.h
PViewOptions.o: PViewOptions.cpp PViewOptions.h ColorTable.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/GmshMessage.h
@@ -138,8 +139,8 @@ OctreePost.o: OctreePost.cpp ../Common/Octree.h \
../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Numeric/Numeric.h \
../Numeric/NumericEmbedded.h ../Common/GmshMessage.h shapeFunctions.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
- ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h ../Geo/SVector3.h \
- ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
+ ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h ../Geo/SVector3.h \
../Numeric/FunctionSpace.h ../Common/GmshMatrix.h
ColorTable.o: ColorTable.cpp ../Common/GmshMessage.h ColorTable.h \
../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/PartitionOptions.h \
--
GitLab