From cc86cd4942d19adccc94b039954dc35702f3cdba Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 15 Apr 2009 14:07:11 +0000
Subject: [PATCH] *** empty log message ***
---
Common/Makefile | 22 +-
Fltk/Makefile | 27 +-
Fltk/menuWindow.cpp | 2 +-
Geo/GEdge.cpp | 2 +-
Geo/GEdge.h | 7 +-
Geo/GEdgeCompound.cpp | 7 +-
Geo/GFaceCompound.cpp | 155 ++--
Geo/GFaceCompound.h | 16 +-
Geo/MElement.cpp | 13 +-
Geo/Makefile | 16 +-
Geo/STensor3.h | 139 ++++
Mesh/BackgroundMesh.cpp | 51 +-
Mesh/Field.cpp | 82 +-
Mesh/Field.h | 6 +-
Mesh/HighOrder.cpp | 150 ++--
Mesh/HighOrder.h | 13 +
Mesh/Makefile | 76 +-
Mesh/gmshSmoothHighOrder.cpp | 1105 ++++++++++++++++++++++++---
Mesh/gmshSmoothHighOrder.h | 61 +-
Mesh/meshGEdge.cpp | 38 +-
Mesh/meshGFace.cpp | 21 +
Mesh/meshGFaceBDS.cpp | 23 +-
Mesh/meshGFaceDelaunayInsertion.cpp | 22 +-
Mesh/qualityMeasures.cpp | 23 +-
Numeric/GmshMatrix.cpp | 60 ++
Numeric/GmshMatrix.h | 27 +
Numeric/Numeric.cpp | 122 +++
Numeric/Numeric.h | 2 +
Numeric/gmshAssembler.h | 7 +
Numeric/gmshElasticity.cpp | 2 +-
Numeric/gmshElasticity.h | 3 +-
Numeric/gmshLaplace.cpp | 4 +
Numeric/gmshLaplace.h | 13 +-
Numeric/gmshLinearSystemGmm.h | 4 +-
Numeric/gmshTermOfFormulation.h | 3 +-
Parser/Makefile | 6 +-
Plugin/Makefile | 37 +-
Post/shapeFunctions.h | 1 +
benchmarks/2d/naca12_2d.geo | 20 +-
39 files changed, 2018 insertions(+), 370 deletions(-)
create mode 100644 Geo/STensor3.h
diff --git a/Common/Makefile b/Common/Makefile
index e615364441..36e3c13d28 100644
--- a/Common/Makefile
+++ b/Common/Makefile
@@ -67,7 +67,9 @@ Gmsh${OBJEXT}: Gmsh.cpp GmshConfig.h GmshDefines.h ../Numeric/GmshPredicates.h \
../Geo/Pair.h ../Geo/GRegion.h ../Geo/GEntity.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h GmshMessage.h OpenFile.h CreateFile.h Options.h \
../Post/ColorTable.h CommandLine.h OS.h ../Mesh/Generator.h \
- ../Mesh/Field.h ../Common/GmshConfig.h ../Post/PView.h Context.h \
+ ../Mesh/Field.h ../Common/GmshConfig.h ../Geo/STensor3.h \
+ ../Geo/SVector3.h ../Numeric/GmshMatrix.h ../Common/GmshMessage.h \
+ ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Post/PView.h Context.h \
../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
../Mesh/meshPartition.h GmshDaemon.h ../Plugin/PluginManager.h
GmshMessage${OBJEXT}: GmshMessage.cpp GmshConfig.h GmshMessage.h Gmsh.h \
@@ -88,14 +90,16 @@ Options${OBJEXT}: Options.cpp GmshConfig.h GmshDefines.h ../Geo/GModel.h \
../Geo/SBoundingBox3d.h GmshMessage.h ../Fltk/Draw.h \
../Mesh/Generator.h Context.h ../Geo/CGNSOptions.h \
../Mesh/meshPartitionOptions.h Options.h ../Post/ColorTable.h \
- DefaultOptions.h ../Mesh/Field.h ../Common/GmshConfig.h ../Post/PView.h \
- ../Mesh/BackgroundMesh.h ../Post/PViewOptions.h ../Post/ColorTable.h \
- ../Post/PViewData.h ../Numeric/GmshMatrix.h ../Common/GmshMessage.h \
- ../Post/adaptiveData.h ../Plugin/PluginManager.h ../Plugin/Plugin.h \
- ../Common/Options.h ../Post/PViewDataList.h ../Post/PViewData.h \
- ../Common/ListUtils.h ../Fltk/GUI.h ../Fltk/Solvers.h \
- ../Fltk/menuWindow.h ../Fltk/popupButton.h ../Fltk/graphicWindow.h \
- ../Fltk/openglWindow.h ../Graphics/drawContext.h ../Fltk/optionWindow.h \
+ DefaultOptions.h ../Mesh/Field.h ../Common/GmshConfig.h \
+ ../Geo/STensor3.h ../Geo/SVector3.h ../Numeric/GmshMatrix.h \
+ ../Common/GmshMessage.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
+ ../Post/PView.h ../Mesh/BackgroundMesh.h ../Post/PViewOptions.h \
+ ../Post/ColorTable.h ../Post/PViewData.h ../Post/adaptiveData.h \
+ ../Plugin/PluginManager.h ../Plugin/Plugin.h ../Common/Options.h \
+ ../Post/PViewDataList.h ../Post/PViewData.h ../Common/ListUtils.h \
+ ../Fltk/GUI.h ../Fltk/Solvers.h ../Fltk/menuWindow.h \
+ ../Fltk/popupButton.h ../Fltk/graphicWindow.h ../Fltk/openglWindow.h \
+ ../Graphics/drawContext.h ../Fltk/optionWindow.h \
../Fltk/spherePositionWidget.h ../Fltk/colorbarWindow.h \
../Fltk/solverWindow.h ../Fltk/manipWindow.h ../Fltk/messageWindow.h \
../Fltk/contextWindow.h ../Fltk/clippingWindow.h
diff --git a/Fltk/Makefile b/Fltk/Makefile
index eedf5f0376..d6c3571841 100644
--- a/Fltk/Makefile
+++ b/Fltk/Makefile
@@ -82,7 +82,9 @@ Main${OBJEXT}: Main.cpp GUI.h menuWindow.h popupButton.h ../Common/Gmsh.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 ../Mesh/Field.h ../Common/GmshConfig.h \
- ../Post/PView.h ../Mesh/BackgroundMesh.h
+ ../Geo/STensor3.h ../Geo/SVector3.h ../Numeric/GmshMatrix.h \
+ ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Post/PView.h \
+ ../Mesh/BackgroundMesh.h
GUI${OBJEXT}: GUI.cpp GUI.h graphicWindow.h openglWindow.h \
../Graphics/drawContext.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
menuWindow.h popupButton.h optionWindow.h spherePositionWidget.h \
@@ -102,10 +104,12 @@ GUI${OBJEXT}: GUI.cpp GUI.h graphicWindow.h openglWindow.h \
../Geo/MVertex.h ../Geo/SVector3.h ../Geo/MFace.h ../Geo/MVertex.h \
../Geo/MEdge.h ../Geo/SVector3.h ../Numeric/FunctionSpace.h \
../Numeric/GmshMatrix.h ../Numeric/Gauss.h ../Post/PView.h Solvers.h \
- ../Mesh/Field.h ../Plugin/Plugin.h ../Common/Options.h \
- ../Post/PViewDataList.h ../Post/PViewData.h ../Common/ListUtils.h \
- ../Plugin/PluginManager.h ../Common/OpenFile.h Win32Icon.h \
- ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
+ ../Mesh/Field.h ../Geo/STensor3.h ../Geo/SVector3.h \
+ ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Plugin/Plugin.h \
+ ../Common/Options.h ../Post/PViewDataList.h ../Post/PViewData.h \
+ ../Common/ListUtils.h ../Plugin/PluginManager.h ../Common/OpenFile.h \
+ Win32Icon.h ../Common/Context.h ../Geo/CGNSOptions.h \
+ ../Mesh/meshPartitionOptions.h
graphicWindow${OBJEXT}: graphicWindow.cpp GUI.h graphicWindow.h openglWindow.h \
../Graphics/drawContext.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
paletteWindow.h mainWindow.h menuWindow.h popupButton.h messageWindow.h \
@@ -154,10 +158,11 @@ menuWindow${OBJEXT}: menuWindow.cpp ../Common/GmshConfig.h \
Solvers.h ../Common/CommandLine.h ../Mesh/Generator.h \
../Mesh/HighOrder.h ../Post/PView.h ../Post/PViewData.h \
../Post/PViewOptions.h ../Post/ColorTable.h ../Mesh/Field.h \
- ../Common/OS.h ../Common/StringUtils.h ../Common/OpenFile.h \
- ../Common/CreateFile.h ../Geo/findLinks.h ../Common/ListUtils.h \
- ../Geo/GeoStringInterface.h ../Common/Context.h ../Geo/CGNSOptions.h \
- ../Mesh/meshPartitionOptions.h
+ ../Geo/STensor3.h ../Geo/SVector3.h ../Numeric/Numeric.h \
+ ../Numeric/GmshMatrix.h ../Common/OS.h ../Common/StringUtils.h \
+ ../Common/OpenFile.h ../Common/CreateFile.h ../Geo/findLinks.h \
+ ../Common/ListUtils.h ../Geo/GeoStringInterface.h ../Common/Context.h \
+ ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
optionWindow${OBJEXT}: optionWindow.cpp ../Common/GmshConfig.h \
../Common/GmshDefines.h ../Common/GmshMessage.h GUI.h optionWindow.h \
spherePositionWidget.h colorbarWindow.h ../Post/ColorTable.h \
@@ -190,7 +195,9 @@ fieldWindow${OBJEXT}: fieldWindow.cpp GUI.h Draw.h fieldWindow.h paletteWindow.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 ../Post/PView.h ../Common/GmshMessage.h \
- ../Mesh/Field.h ../Common/GmshConfig.h ../Geo/GeoStringInterface.h \
+ ../Mesh/Field.h ../Common/GmshConfig.h ../Geo/STensor3.h \
+ ../Geo/SVector3.h ../Numeric/GmshMatrix.h ../Numeric/Numeric.h \
+ ../Numeric/GmshMatrix.h ../Geo/GeoStringInterface.h \
../Common/ListUtils.h ../Common/Options.h ../Post/ColorTable.h \
../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
pluginWindow${OBJEXT}: pluginWindow.cpp GUI.h Draw.h pluginWindow.h \
diff --git a/Fltk/menuWindow.cpp b/Fltk/menuWindow.cpp
index 9f744901e6..1afbfd1dc6 100644
--- a/Fltk/menuWindow.cpp
+++ b/Fltk/menuWindow.cpp
@@ -2410,7 +2410,7 @@ contextItem menu_mesh[] = {
#if defined(HAVE_FOURIER_MODEL)
{"Reparameterize", (Fl_Callback *)mesh_parameterize_cb} ,
#endif
- // {"Reclassify", (Fl_Callback *)mesh_classify_cb} ,
+ {"Reclassify", (Fl_Callback *)mesh_classify_cb} ,
{"Save", (Fl_Callback *)mesh_save_cb} ,
{""}
};
diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp
index 85309dda0a..e1379d1a36 100644
--- a/Geo/GEdge.cpp
+++ b/Geo/GEdge.cpp
@@ -21,7 +21,7 @@
#endif
GEdge::GEdge(GModel *model, int tag, GVertex *_v0, GVertex *_v1)
- : GEntity(model, tag), _tooSmall(false), v0(_v0), v1(_v1)
+ : GEntity(model, tag), _tooSmall(false), v0(_v0), v1(_v1), compound(0)
{
if(v0) v0->addEdge(this);
if(v1 && v1 != v0) v1->addEdge(this);
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index a4489807ac..3cdad01497 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -16,6 +16,7 @@
class MElement;
class MLine;
class ExtrudeParams;
+class GEdgeCompound;
#include <set>
@@ -27,8 +28,8 @@ class GEdge : public GEntity {
protected:
GVertex *v0, *v1;
+ GEdgeCompound *compound; // this model edge belongs to a compound
std::list<GFace *> l_faces;
-
std::set<GFace *> bl_faces;
public:
@@ -153,6 +154,10 @@ class GEdge : public GEntity {
virtual bool XYZToU(const SVector3& Q,double& t,const double relax=0.5) const;
+ // compound
+ void setCompound (GEdgeCompound *gec) {compound = gec;}
+ GEdgeCompound *getCompound () const {return compound;}
+
struct {
char Method;
double coeffTransfinite;
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index aba627092c..f2ad5d907a 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -23,13 +23,16 @@ void GEdgeCompound::orderEdges()
{
std::list<GEdge*> edges ;
std::vector<GEdge*> _c ;
- for (int i=0;i<_compound.size();i++)edges.push_back(_compound[i]);
+ for (int i=0;i<_compound.size();i++){
+ _compound[i]->setCompound(this);
+ edges.push_back(_compound[i]);
+ }
_c.push_back(*(edges.begin()));
edges.erase(edges.begin());
_orientation.push_back(true);
GVertex *first = _c[0]->getBeginVertex();
GVertex *last = _c[0]->getEndVertex();
-
+
while (first != last){
if (edges.empty())break;
for (std::list<GEdge*>::iterator it = edges.begin() ; it != edges.end() ; ++it){
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 1a0542bb1f..6268d4a1b2 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -45,6 +45,29 @@ class gmshGradientBasedDiffusivity : public gmshFunction<double>
}
};
+class gmshDistanceBasedDiffusivity : public gmshFunction<double>
+{
+ private:
+ MElement *_current;
+ mutable std::map<MVertex*, SPoint3> _coordinates;
+ public:
+ gmshDistanceBasedDiffusivity (std::map<MVertex*,SPoint3> &coordinates)
+ : _current (0),_coordinates(coordinates){}
+ void setCurrent (MElement *current){ _current = current; }
+ virtual double operator () (double x, double y, double z) const
+ {
+ double xyz[3] = {x, y, z}, uvw[3];
+ _current->xyz2uvw(xyz, uvw);
+ double value[256];
+ for (int i = 0; i < _current->getNumVertices(); i++){
+ const SPoint3 p = _coordinates[_current->getVertex(i)];
+ value[i] = p[2];
+ }
+ double val = _current->interpolate(value, uvw[0], uvw[1], uvw[2]);
+ return 1.0;//exp(15*val);
+ }
+};
+
static void fixEdgeToValue(GEdge *ed, double value, gmshAssembler<double> &myAssembler)
{
myAssembler.fixVertex(ed->getBeginVertex()->mesh_vertices[0], 0, 1, value);
@@ -58,10 +81,9 @@ void GFaceCompound::parametrize() const
{
if (!oct){
coordinates.clear();
- parametrize(0, 0);
- parametrize(1, 0);
- // parametrize(0,1);
- // parametrize(1,1);
+ parametrize(ITERD);
+ parametrize(ITERU);
+ parametrize(ITERV);
computeNormals();
buildOct();
}
@@ -69,17 +91,12 @@ void GFaceCompound::parametrize() const
void GFaceCompound::getBoundingEdges()
{
- if (_U0.size() && !_V1.size()){
+ if (_U0.size()){
std::list<GEdge*> :: const_iterator it = _U0.begin();
for ( ; it != _U0.end() ; ++it){
l_edges.push_back(*it);
(*it)->addFace(this);
}
- it = _U1.begin();
- for ( ; it != _U1.end() ; ++it){
- l_edges.push_back(*it);
- (*it)->addFace(this);
- }
it = _V0.begin();
for ( ; it != _V0.end() ; ++it){
l_edges.push_back(*it);
@@ -112,6 +129,8 @@ void GFaceCompound::getBoundingEdges()
l_edges.push_back(*itf);
(*itf)->addFace(this);
}
+
+ _U0 = l_edges;
}
GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
@@ -121,9 +140,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
{
getBoundingEdges();
if (!_U0.size()) _type = UNITCIRCLE;
- else if (!_V0.size()) _type = UNITCIRCLE;
- else if (!_U1.size()) _type = CYLINDER;
- else if (!_V1.size()) _type = BIFURCATION;
+ else if (!_V1.size()) _type = UNITCIRCLE;
else _type = SQUARE;
}
@@ -196,26 +213,29 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
return true;
}
-void GFaceCompound::parametrize(bool _isU, int ITER) const
+void GFaceCompound::parametrize(iterationStep step) const
{
- Msg::Info("Parametrizing Surface %d coordinate %d (ITER %d)", tag(), _isU, ITER);
+ Msg::Info("Parametrizing Surface %d coordinate (ITER %d)", tag(), step);
#ifdef HAVE_GMM
gmshLinearSystemGmm<double> lsys;
- lsys.setPrec(1.e-10);
+ lsys.setPrec(1.e-15);
if(Msg::GetVerbosity() == 99) lsys.setNoisy(2);
#else
gmshLinearSystemFull<double> lsys;
#endif
gmshAssembler<double> myAssembler(&lsys);
- gmshGradientBasedDiffusivity diffusivity(coordinates);
- if (ITER > 0) diffusivity.setComponent(_isU);
-
+
+ gmshDistanceBasedDiffusivity diffusivity(coordinates);
+ gmshFunction<double> ONE(1.0);
+
if (_type == UNITCIRCLE){
// maps the boundary onto a circle
std::vector<MVertex*> ordered;
+ std::vector<MVertex*> ordered1;
std::vector<double> coords;
- bool success = orderVertices(l_edges, ordered, coords);
+ std::vector<double> coords1;
+ bool success = orderVertices(_U0, ordered, coords);
if (!success){
Msg::Error("Could not order vertices on boundary");
return;
@@ -223,25 +243,24 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
for (unsigned int i = 0; i < ordered.size(); i++){
MVertex *v = ordered[i];
const double theta = 2 * M_PI * coords[i];
- if (_isU) myAssembler.fixVertex(v, 0, 1, cos(theta));
- else myAssembler.fixVertex(v, 0, 1, sin(theta));
+ if (step == ITERU) myAssembler.fixVertex(v, 0, 1, cos(theta));
+ else if (step == ITERV) myAssembler.fixVertex(v, 0, 1, sin(theta));
+ else if (step == ITERD) myAssembler.fixVertex(v, 0, 1, 1.0);
}
- }
- else if (_type == CYLINDER){
- // maps the boundary onto an annulus
- std::vector<MVertex*> ordered;
- std::vector<double> coords;
- bool success = orderVertices (_U0, ordered, coords);
- if (!success)throw;
- for (unsigned int i = 0; i < ordered.size(); i++){
- MVertex *v = ordered[i];
- const double theta = 2 * M_PI * coords[i];
- if (_isU) myAssembler.fixVertex(v, 0, 1, cos(theta));
- else myAssembler.fixVertex(v, 0, 1, sin(theta));
+ if (step == ITERD && _V0.size()) {
+ success = orderVertices(_V0, ordered1, coords1);
+ if (!success){
+ Msg::Error("Could not order vertices on boundary");
+ return;
+ }
+ for (unsigned int i = 0; i < ordered1.size(); i++){
+ MVertex *v = ordered1[i];
+ myAssembler.fixVertex(v, 0, 1, 0.0);
+ }
}
}
- else{
- if (_isU){
+ else if (_type == SQUARE){
+ if (step == ITERU){
std::list<GEdge*> :: const_iterator it = _U0.begin();
for ( ; it != _U0.end(); ++it){
fixEdgeToValue(*it, 0.0, myAssembler);
@@ -251,7 +270,7 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
fixEdgeToValue(*it, 1.0, myAssembler);
}
}
- else{
+ else if (step == ITERV){
std::list<GEdge*> :: const_iterator it = _V0.begin();
for ( ; it != _V0.end(); ++it){
fixEdgeToValue (*it, 0.0, myAssembler);
@@ -262,7 +281,8 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
}
}
}
-
+ else throw;
+
std::list<GFace*>::const_iterator it = _compound.begin();
for ( ; it != _compound.end(); ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -276,15 +296,28 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
Msg::Debug("Creating term %d dofs numbered %d fixed",
myAssembler.sizeOfR(), myAssembler.sizeOfF());
- gmshLaplaceTerm laplace(model(), &diffusivity, 1);
- it = _compound.begin();
- for ( ; it != _compound.end() ; ++it){
- for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- MTriangle *t = (*it)->triangles[i];
- diffusivity.setCurrent(t);
- laplace.addToMatrix(myAssembler, t);
- }
- }
+ if (step == ITERD){
+ gmshLaplaceTerm laplace(model(), &ONE, 1);
+ it = _compound.begin();
+ for ( ; it != _compound.end() ; ++it){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ MTriangle *t = (*it)->triangles[i];
+ diffusivity.setCurrent(t);
+ laplace.addToMatrix(myAssembler, t);
+ }
+ }
+ }
+ else{
+ gmshLaplaceTerm laplace(model(), &diffusivity, 1);
+ it = _compound.begin();
+ for ( ; it != _compound.end() ; ++it){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
+ MTriangle *t = (*it)->triangles[i];
+ diffusivity.setCurrent(t);
+ laplace.addToMatrix(myAssembler, t);
+ }
+ }
+ }
Msg::Debug("Assembly done");
lsys.systemSolve();
@@ -303,14 +336,14 @@ void GFaceCompound::parametrize(bool _isU, int ITER) const
for (std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
double value = myAssembler.getDofValue(v,0,1);
- std::map<MVertex*,SPoint2>::iterator itf = coordinates.find(v);
+ std::map<MVertex*,SPoint3>::iterator itf = coordinates.find(v);
if (itf == coordinates.end()){
- SPoint2 p (_isU ? value:0.0,_isU ? 0.0:value);
+ SPoint3 p(0,0,0);
+ p[step] = value;
coordinates[v] = p;
}
else{
- if(_isU) itf->second[0]= value;
- else itf->second[1]= value;
+ itf->second[step]= value;
}
}
}
@@ -453,9 +486,9 @@ static int GFaceCompoundInEle(void *a, double*c)
GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a;
double M[2][2],R[2],X[2];
const double eps = 1.e-8;
- const SPoint2 p0 = t->p1;
- const SPoint2 p1 = t->p2;
- const SPoint2 p2 = t->p3;
+ const SPoint3 p0 = t->p1;
+ const SPoint3 p1 = t->p2;
+ const SPoint3 p2 = t->p3;
M[0][0] = p1.x() - p0.x();
M[0][1] = p2.x() - p0.x();
M[1][0] = p1.y() - p0.y();
@@ -480,9 +513,9 @@ void GFaceCompound::getTriangle(double u, double v,
if (!(*lt)) return;
double M[2][2],X[2],R[2];
- const SPoint2 p0 = (*lt)->p1;
- const SPoint2 p1 = (*lt)->p2;
- const SPoint2 p2 = (*lt)->p3;
+ const SPoint3 p0 = (*lt)->p1;
+ const SPoint3 p1 = (*lt)->p2;
+ const SPoint3 p2 = (*lt)->p3;
M[0][0] = p1.x() - p0.x();
M[0][1] = p2.x() - p0.x();
M[1][0] = p1.y() - p0.y();
@@ -504,7 +537,7 @@ void GFaceCompound::buildOct() const
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
for (int j=0;j<3;j++){
- std::map<MVertex*,SPoint2>::const_iterator itj =
+ std::map<MVertex*,SPoint3>::const_iterator itj =
coordinates.find(t->getVertex(j));
bb += SPoint3(itj->second.x(),itj->second.y(),0.0);
}
@@ -539,11 +572,11 @@ void GFaceCompound::buildOct() const
for ( ; it != _compound.end() ; ++it){
for (unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
- std::map<MVertex*,SPoint2>::const_iterator it0 =
+ std::map<MVertex*,SPoint3>::const_iterator it0 =
coordinates.find(t->getVertex(0));
- std::map<MVertex*,SPoint2>::const_iterator it1 =
+ std::map<MVertex*,SPoint3>::const_iterator it1 =
coordinates.find(t->getVertex(1));
- std::map<MVertex*,SPoint2>::const_iterator it2 =
+ std::map<MVertex*,SPoint3>::const_iterator it2 =
coordinates.find(t->getVertex(2));
_gfct[count].p1 = it0->second;
_gfct[count].p2 = it1->second;
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 5ceae9b7d3..669fe7a202 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -28,7 +28,7 @@ generator of gmsh!
class GFaceCompoundTriangle {
public:
- SPoint2 p1, p2, p3;
+ SPoint3 p1, p2, p3;
SPoint2 gfp1, gfp2, gfp3;
SPoint3 v1, v2, v3;
GFace *gf;
@@ -39,23 +39,25 @@ public:
class Octree;
class GFaceCompound : public GFace {
+ public:
+ typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
protected:
std::list<GFace*> _compound;
std::list<GEdge*> _U0, _U1, _V0, _V1;
mutable GFaceCompoundTriangle *_gfct;
mutable Octree *oct;
- mutable std::map<MVertex*,SPoint2> coordinates;
+ mutable std::map<MVertex*,SPoint3> coordinates;
mutable std::map<MVertex*,SVector3> _normals;
void buildOct() const ;
- void parametrize(bool,int) const ;
void parametrize() const ;
+ void parametrize(iterationStep) const ;
void computeNormals () const;
void getBoundingEdges();
void getTriangle(double u, double v, GFaceCompoundTriangle **lt,
double &_u, double &_v) const;
virtual double curvature(MTriangle *t) const;
public:
- typedef enum {UNITCIRCLE, CYLINDER, BIFURCATION, SQUARE} typeOfIsomorphism;
+ typedef enum {UNITCIRCLE, SQUARE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag,
std::list<GFace*> &compound,
std::list<GEdge*> &U0,
@@ -69,7 +71,11 @@ public:
virtual GEntity::GeomType geomType() const { return CompoundSurface; }
ModelType getNativeType() const { return GmshModel; }
void * getNativePtr() const { return 0; }
- SPoint2 getCoordinates(MVertex *v) const { parametrize() ; return coordinates[v]; }
+ SPoint2 getCoordinates(MVertex *v) const {
+ parametrize() ;
+ std::map<MVertex*,SPoint3>::iterator it = coordinates.find(v);
+ return SPoint2(it->second.x(),it->second.y());
+ }
virtual bool buildRepresentationCross(){ return false; }
virtual double curvature(const SPoint2 ¶m) const;
private:
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index d5689b65d9..fa48c1dab8 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -144,9 +144,16 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
switch (ele->getDim()) {
+ case 0:
+ {
+ jac[0][0] = jac[1][1] = jac[2][2] = 1.0;
+ jac[0][1] = jac[1][0] = jac[2][0] = 0.0;
+ jac[0][2] = jac[1][2] = jac[2][1] = 0.0;
+ break;
+ }
case 1:
{
- dJ = sqrt(SQU(jac[0][0]) + SQU(jac[0][1]) + SQU(jac[0][2]));
+ dJ = sqrt(SQU(jac[0][0]) + SQU(jac[1][0]) + SQU(jac[2][0]));
// regularize matrix
double a[3], b[3], c[3];
@@ -161,8 +168,8 @@ static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
b[0] = 0.; b[1] = a[2]; b[2] = -a[1];
}
prodve(a, b, c);
- jac[1][0] = b[0]; jac[1][1] = b[1]; jac[1][2] = b[2];
- jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
+ jac[0][1] = b[0]; jac[1][1] = b[1]; jac[2][1] = b[2];
+ jac[0][2] = c[0]; jac[1][2] = c[1]; jac[2][2] = c[2];
break;
}
case 2:
diff --git a/Geo/Makefile b/Geo/Makefile
index 28fc4a8679..acc4f2e589 100644
--- a/Geo/Makefile
+++ b/Geo/Makefile
@@ -236,9 +236,9 @@ GModel${OBJEXT}: GModel.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
MHexahedron.h MPrism.h MPyramid.h discreteRegion.h discreteFace.h \
discreteEdge.h discreteVertex.h gmshSurface.h ../Numeric/Numeric.h \
../Numeric/GmshMatrix.h ../Common/Octree.h ../Common/OctreeInternals.h \
- ../Common/SmoothData.h ../Mesh/Field.h ../Post/PView.h ../Geo/SPoint3.h \
- ../Mesh/Generator.h ../Common/Context.h ../Geo/CGNSOptions.h \
- ../Mesh/meshPartitionOptions.h
+ ../Common/SmoothData.h ../Mesh/Field.h ../Geo/STensor3.h \
+ ../Geo/SVector3.h ../Post/PView.h ../Geo/SPoint3.h ../Mesh/Generator.h \
+ ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
GModelIO_Geo${OBJEXT}: GModelIO_Geo.cpp ../Common/GmshConfig.h \
../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 \
@@ -247,8 +247,9 @@ GModelIO_Geo${OBJEXT}: GModelIO_Geo.cpp ../Common/GmshConfig.h \
../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
../Common/ListUtils.h ExtrudeParams.h ../Common/SmoothData.h \
../Common/OpenFile.h gmshVertex.h gmshFace.h GFaceCompound.h \
- GEdgeCompound.h gmshEdge.h gmshRegion.h ../Mesh/Field.h ../Post/PView.h \
- ../Geo/SPoint3.h ../Parser/Parser.h
+ GEdgeCompound.h gmshEdge.h gmshRegion.h ../Mesh/Field.h \
+ ../Geo/STensor3.h ../Geo/SVector3.h ../Post/PView.h ../Geo/SPoint3.h \
+ ../Parser/Parser.h
GModelIO_Mesh${OBJEXT}: GModelIO_Mesh.cpp 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/GmshDefines.h MPoint.h \
@@ -292,8 +293,9 @@ Geo${OBJEXT}: Geo.cpp ../Common/GmshMessage.h ../Numeric/Numeric.h \
../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
ExtrudeParams.h ../Common/SmoothData.h GModel.h GVertex.h GEntity.h \
GPoint.h GEdge.h GFace.h GEdgeLoop.h GRegion.h GeoInterpolation.h \
- ../Mesh/Field.h ../Post/PView.h ../Geo/SPoint3.h ../Common/Context.h \
- ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
+ ../Mesh/Field.h ../Geo/STensor3.h ../Geo/SVector3.h ../Post/PView.h \
+ ../Geo/SPoint3.h ../Common/Context.h ../Geo/CGNSOptions.h \
+ ../Mesh/meshPartitionOptions.h
GeoStringInterface${OBJEXT}: GeoStringInterface.cpp ../Common/GmshConfig.h \
../Common/GmshMessage.h ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
../Common/StringUtils.h Geo.h ../Common/GmshDefines.h gmshSurface.h \
diff --git a/Geo/STensor3.h b/Geo/STensor3.h
new file mode 100644
index 0000000000..a4659ffed5
--- /dev/null
+++ b/Geo/STensor3.h
@@ -0,0 +1,139 @@
+// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+
+#ifndef _STENSOR3_H_
+#define _STENSOR3_H_
+
+#include "SVector3.h"
+#include "GmshMatrix.h"
+#include "Numeric.h"
+// concrete class for symmetric positive definite 3x3 matrix
+class SMetric3 {
+ protected:
+ // lower diagonal storage
+ // 00 10 11 20 21 22
+ double _val[6];
+ public:
+ inline int getIndex(int i, int j) const{
+ static int _index[3][3] = {{0,1,3},{1,2,4},{3,4,5}};
+ return _index[i][j];
+ }
+ void getMat (gmshMatrix<double> & mat) const{
+ for (int i=0;i<3;i++)
+ for (int j=0;j<3;j++)
+ mat(i,j) = _val[getIndex(i,j)];
+ }
+ void setMat (const gmshMatrix<double> & mat){
+ for (int i=0;i<3;i++)
+ for (int j=0;j<3;j++)
+ _val[getIndex(i,j)] = mat(i,j);
+ }
+ // default constructor, identity
+ SMetric3(const double v = 1.0){
+ _val[0] = _val[2] = _val[5] = v;
+ _val[1] = _val[3] = _val[4] = 0.0;
+ }
+ SMetric3(const double l1,
+ const double l2,
+ const double l3,
+ const SVector3 &t1,
+ const SVector3 &t2,
+ const SVector3 &t3){
+ SVector3 t1b = t1;
+ SVector3 t2b = t2;
+ SVector3 t3b = t3;
+ t1b[0] *= l1; t2b[0] *= l1; t3b[0] *= l1;
+ t1b[1] *= l2; t2b[1] *= l2; t3b[1] *= l2;
+ t1b[2] *= l3; t2b[2] *= l3; t3b[2] *= l3;
+ _val[0] = dot (t1b,t1);
+ _val[1] = dot (t2b,t1);
+ _val[2] = dot (t2b,t2);
+ _val[3] = dot (t3b,t1);
+ _val[4] = dot (t3b,t2);
+ _val[5] = dot (t3b,t3);
+ }
+ inline double &operator()(int i, int j){
+ return _val[getIndex(i,j)];
+ }
+ inline double operator()(int i, int j) const{
+ return _val[getIndex(i,j)];
+ }
+ SMetric3 invert () const {
+ gmshMatrix<double> m(3,3);
+ getMat(m);
+ m.invertInPlace();
+ SMetric3 ithis;
+ ithis.setMat(m);
+ return ithis;
+ }
+ SMetric3 operator + (const SMetric3 &other) const {
+ SMetric3 res(*this);
+ for (int i=0;i<6;i++)res._val[i]+=other._val[i];
+ return res;
+ }
+ SMetric3& operator += (const SMetric3 &other) {
+ for (int i=0;i<6;i++)_val[i]+=other._val[i];
+ return *this;
+ }
+ SMetric3& operator *= (const double &other) {
+ for (int i=0;i<6;i++)_val[i]*=other;
+ return *this;
+ }
+ SMetric3& operator *= (const SMetric3 &other) {
+ gmshMatrix<double> m1(3,3),m2(3,3),m3(3,3);
+ getMat(m1);
+ other.getMat(m2);
+ m1.mult(m2,m3);
+ setMat(m3);
+ return *this;
+ }
+ void eig (gmshMatrix<double> &V, gmshVector<double> &S) {
+ gmshMatrix<double> me(3,3),right(3,3);
+ gmshVector<double> im(3);
+ getMat(me);
+ me.eig(V,S,im,right);
+ }
+};
+
+// scalar product with respect to the metric
+inline double dot(const SVector3 &a, const SMetric3 &m, const SVector3 &b)
+{ return
+ b.x() * ( m(0,0) * a.x() + m(0,1) * a.y() + m(0,2) * a.z() ) +
+ b.y() * ( m(1,0) * a.x() + m(1,1) * a.y() + m(1,2) * a.z() ) +
+ b.z() * ( m(2,0) * a.x() + m(2,1) * a.y() + m(2,2) * a.z() ) ;}
+
+// compute the largest inscribed ellipsoid...
+inline SMetric3 intersection (const SMetric3 &m1, const SMetric3 &m2)
+{
+ SMetric3 im1 = m1.invert();
+ gmshMatrix<double> V(3,3);
+ gmshVector<double> S(3);
+ im1 *= m2;
+ im1.eig(V,S);
+ SVector3 v0(V(0,0),V(1,0),V(2,0));
+ SVector3 v1(V(0,1),V(1,1),V(2,1));
+ SVector3 v2(V(0,2),V(1,2),V(2,2));
+ double l0 = std::max(dot(v0,m1,v0),dot(v0,m2,v0));
+ double l1 = std::max(dot(v1,m1,v1),dot(v1,m2,v1));
+ double l2 = std::max(dot(v2,m1,v2),dot(v2,m2,v2));
+
+ return SMetric3(l0,l1,l2,v0,v1,v2);
+}
+
+// (1-t) * m1 + t * m2
+inline SMetric3 interpolation (const SMetric3 &m1,
+ const SMetric3 &m2,
+ const double t)
+{
+ SMetric3 im1 = m1.invert();
+ SMetric3 im2 = m2.invert();
+ im1 *= (1.-t);
+ im2 *= t;
+ im1 += im2;
+ return im1.invert();
+}
+
+
+#endif
diff --git a/Mesh/BackgroundMesh.cpp b/Mesh/BackgroundMesh.cpp
index 6c05841e43..606605cc32 100644
--- a/Mesh/BackgroundMesh.cpp
+++ b/Mesh/BackgroundMesh.cpp
@@ -135,7 +135,8 @@ static double LC_MVertex_PNTS(GEntity *ge, double U, double V)
}
// This is the only function that is used by the meshers
-double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double Z)
+double BGM_MeshSize(GEntity *ge, double U, double V,
+ double X, double Y, double Z)
{
// default lc (mesh size == size of the model)
double l1 = CTX::instance()->lc;
@@ -172,6 +173,54 @@ double BGM_MeshSize(GEntity *ge, double U, double V, double X, double Y, double
return lc * CTX::instance()->mesh.lcFactor;
}
+
+// anisotropic version
+SMetric3 BGM_MeshMetric(GEntity *ge,
+ double U, double V,
+ double X, double Y, double Z)
+{
+ // default lc (mesh size == size of the model)
+ SMetric3 l1(CTX::instance()->lc);
+
+ // lc from points
+ SMetric3 l2(MAX_LC);
+ if(CTX::instance()->mesh.lcFromPoints && ge->dim() < 2)
+ l2 = SMetric3(LC_MVertex_PNTS(ge, U, V));
+
+ // lc from curvature
+ SMetric3 l3 (MAX_LC);
+ if(CTX::instance()->mesh.lcFromCurvature && ge->dim() < 3)
+ l3 = SMetric3(LC_MVertex_CURV(ge, U, V));
+
+ // lc from fields
+ SMetric3 l4 (MAX_LC);
+ FieldManager *fields = GModel::current()->getFields();
+ if(fields->background_field > 0){
+ Field *f = fields->get(fields->background_field);
+ if(f){
+ // if (!f->isotropic())
+ // (*f)(X, Y, Z, ge,l4);
+ // else
+ l4 = SMetric3((*f)(X, Y, Z, ge));
+ }
+ }
+
+ // take the minimum, then constrain by lcMin and lcMax
+ // double lc = std::min(std::min(std::min(l1, l2), l3), l4);
+
+ // lc = std::max(lc, CTX::instance()->mesh.lcMin);
+ // lc = std::min(lc, CTX::instance()->mesh.lcMax);
+
+ // if(lc <= 0.){
+ /// Msg::Error("Wrong characteristic length lc = %g (lcmin = %g, lcmax = %g)",
+ // lc, CTX::instance()->mesh.lcMin, CTX::instance()->mesh.lcMax);
+ // lc = l1;
+ // }
+
+// return lc * CTX::instance()->mesh.lcFactor;
+}
+
+
bool Extend1dMeshIn2dSurfaces()
{
return CTX::instance()->mesh.lcExtendFromBoundary ? true : false;
diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp
index ffcbf5f828..40abe3b406 100644
--- a/Mesh/Field.cpp
+++ b/Mesh/Field.cpp
@@ -477,15 +477,16 @@ class BoxField : public Field
class ThresholdField : public Field
{
+ protected :
int iField;
double dmin, dmax, lcmin, lcmax;
bool sigmoid, stopAtDistMax;
public:
- const char *getName()
+ virtual const char *getName()
{
return "Threshold";
}
- std::string getDescription()
+ virtual std::string getDescription()
{
return "F = LCMin if Field[IField] <= DistMin,\n"
"F = LCMax if Field[IField] >= DistMax,\n"
@@ -538,6 +539,80 @@ class ThresholdField : public Field
}
};
+class BoundaryLayerField : public ThresholdField {
+public:
+ BoundaryLayerField()
+ { }
+ virtual bool isotropic () const {return false;}
+ virtual const char *getName()
+ {
+ return "BoundaryLayer";
+ }
+ virtual std::string getDescription()
+ {
+ return "F = LCMin if Field[IField] <= DistMin,\n"
+ "F = LCMax if Field[IField] >= DistMax,\n"
+ "F = interpolation between LcMin and LcMax if DistMin < Field[IField] < DistMax";
+ }
+ virtual void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0)
+ {
+ Field *field = GModel::current()->getFields()->get(iField);
+ if(!field) {
+ metr(0,0) = 1/(MAX_LC*MAX_LC);
+ metr(1,1) = 1/(MAX_LC*MAX_LC);
+ metr(2,2) = 1/(MAX_LC*MAX_LC);
+ metr(0,1) = metr(0,2) = metr(1,2) = 0;
+ return;
+ }
+ double dist = (*field) (x, y, z);
+
+ double r = (dist - dmin) / (dmax - dmin);
+ r = std::max(std::min(r, 1.), 0.);
+ double lc;
+ if(stopAtDistMax && r >= 1.){
+ lc = MAX_LC;
+ }
+ else if(sigmoid){
+ double s = exp(12. * r - 6.) / (1. + exp(12. * r - 6.));
+ lc = lcmin * (1. - s) + lcmax * s;
+ }
+ else{ // linear
+ lc = lcmin * (1 - r) + lcmax * r;
+ }
+
+ double delta = std::min(CTX::instance()->lc / 1e4, dist);
+ double gx =
+ ((*field) (x + delta / 2, y, z) -
+ (*field) (x - delta / 2, y, z)) / delta;
+ double gy =
+ ((*field) (x, y + delta / 2, z) -
+ (*field) (x, y - delta / 2, z)) / delta;
+ double gz =
+ ((*field) (x, y, z + delta / 2) -
+ (*field) (x, y, z - delta / 2)) / delta;
+
+ SVector3 g(gx,gy,gz);
+ g.normalize();
+ SVector3 t1,t2;
+
+ if (fabs(gx) < fabs(gy) && fabs(gx) < fabs(gz))
+ t1 = SVector3(1,0,0);
+ else if (fabs(gy) < fabs(gx) && fabs(gy) < fabs(gz))
+ t1 = SVector3(0,1,0);
+ else
+ t1 = SVector3(0,0,1);
+
+ t2 = crossprod(g,t1);
+ t2.normalize();
+ t1 = crossprod(t2,g);
+
+ metr = SMetric3(1./(lc*lc),
+ 1/(lcmax*lcmax),
+ 1/(lcmax*lcmax),
+ g,t1,t2);
+ }
+};
+
class GradientField : public Field
{
int iField, kind;
@@ -600,6 +675,7 @@ class GradientField : public Field
}
};
+
class CurvatureField : public Field
{
int iField;
@@ -1236,6 +1312,7 @@ class AttractorField : public Field
return sqrt(dist[0]);
}
};
+
#endif
template<class F> class FieldFactoryT : public FieldFactory {
@@ -1252,6 +1329,7 @@ FieldManager::FieldManager()
{
map_type_name["Structured"] = new FieldFactoryT<StructuredField>();
map_type_name["Threshold"] = new FieldFactoryT<ThresholdField>();
+ map_type_name["BoundaryLayer"] = new FieldFactoryT<BoundaryLayerField>();
map_type_name["Box"] = new FieldFactoryT<BoxField>();
map_type_name["LonLat"] = new FieldFactoryT<LonLatField>();
#if !defined(HAVE_NO_POST)
diff --git a/Mesh/Field.h b/Mesh/Field.h
index a35e770137..d4ca54e45b 100644
--- a/Mesh/Field.h
+++ b/Mesh/Field.h
@@ -10,6 +10,7 @@
#include <map>
#include <list>
#include "GmshConfig.h"
+#include "STensor3.h"
#if !defined(HAVE_NO_POST)
#include "PView.h"
@@ -62,7 +63,10 @@ class Field {
public:
int id;
std::map<std::string, FieldOption *> options;
- virtual double operator() (double x, double y, double z, GEntity *ge=0) = 0;
+ virtual double operator() (double x, double y, double z, GEntity *ge=0) = 0;
+ // start of the anisotropic field implementation
+ virtual void operator() (double x, double y, double z, SMetric3 &, GEntity *ge=0){throw;}
+ virtual bool isotropic () const {return true;}
virtual ~Field() {}
bool update_needed;
Field();
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 954b6b0c4f..06e100b7e1 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -288,10 +288,10 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
GPoint pc = ge->point(US[j + 1]);
v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge, US[j + 1]);
- if (displ2D){
+ if (displ2D || displ3D){
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()));
+ if (displ2D)displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ if (displ3D)displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
}
}
temp.push_back(v);
@@ -309,9 +309,10 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
}
static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
- edgeContainer &edgeVertices, bool linear,
- int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
- gmshHighOrderSmoother *displ3D = 0)
+ edgeContainer &edgeVertices, bool linear,
+ int nPts = 1,
+ gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0)
{
if(gf->geomType() == GEntity::DiscreteSurface ||
gf->geomType() == GEntity::BoundaryLayerSurface)
@@ -349,9 +350,10 @@ static void getEdgeVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &ve,
else{
GPoint pc = gf->point(US[j + 1], VS[j + 1]);
v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, US[j + 1], VS[j + 1]);
- if (displ3D){
+ if (displ2D || displ3D){
SPoint3 pc2 = edge.interpolate(t);
- displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ if (displ3D) displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ if (displ2D) displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
}
}
temp.push_back(v);
@@ -477,9 +479,10 @@ static void getFaceVertices(GFace *gf, MElement *incomplete, MElement *ele,
else{
v = new MVertex(X, Y, Z, gf);
}
- if(displ3D){
+ if(displ3D || displ2D){
SPoint3 pc2 = face.interpolate(t1, t2);
- displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ if(displ3D)displ3D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
+ if(displ2D)displ2D->add(v, SVector3(pc2.x(), pc2.y(), pc2.z()));
}
}
// should be expensive -> induces a new search each time
@@ -698,6 +701,47 @@ static void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
ge->deleteVertexArrays();
}
+MTriangle* setHighOrder(MTriangle *t,
+ GFace *gf,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ bool linear,
+ bool incomplete,
+ int nPts,
+ gmshHighOrderSmoother *displ2D,
+ gmshHighOrderSmoother *displ3D){
+
+ std::vector<MVertex*> ve, vf;
+ getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts, displ2D, displ3D);
+ if(nPts == 1){
+ return new MTriangle6(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ ve[0], ve[1], ve[2]);
+ }
+ else{
+ if(incomplete){
+ return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ ve, nPts + 1);
+ }
+ else{
+ // MTriangleN incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ // ve, nPts + 1);
+ if (displ2D && gf->geomType() == GEntity::Plane){
+ MTriangle incpl(t->getVertex(0), t->getVertex(1), t->getVertex(2));
+ getFaceVertices(gf, &incpl, t, vf, faceVertices, linear, nPts, displ2D, displ3D);
+ }
+ 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);
+ }
+ ve.insert(ve.end(), vf.begin(), vf.end());
+ return new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
+ ve, nPts + 1);
+ }
+ }
+}
+
+
static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
faceContainer &faceVertices, bool linear, bool incomplete,
int nPts = 1, gmshHighOrderSmoother *displ2D = 0,
@@ -706,29 +750,9 @@ static void setHighOrder(GFace *gf, edgeContainer &edgeVertices,
std::vector<MTriangle*> triangles2;
for(unsigned int i = 0; i < gf->triangles.size(); i++){
MTriangle *t = gf->triangles[i];
- std::vector<MVertex*> ve, vf;
- getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts, displ2D, displ3D);
- if(nPts == 1){
- triangles2.push_back
- (new MTriangle6(t->getVertex(0), t->getVertex(1), t->getVertex(2),
- ve[0], ve[1], ve[2]));
- }
- else{
- if(incomplete){
- triangles2.push_back
- (new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
- ve, nPts + 1));
- }
- 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);
- ve.insert(ve.end(), vf.begin(), vf.end());
- triangles2.push_back
- (new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
- ve, nPts + 1));
- }
- }
+ MTriangle *tNew = setHighOrder(t,gf,edgeVertices,faceVertices, linear, incomplete,
+ nPts,displ2D,displ3D);
+ triangles2.push_back(tNew);
delete t;
}
gf->triangles = triangles2;
@@ -932,29 +956,63 @@ void SetOrder1(GModel *m)
removeHighOrderVertices(*it);
}
-static void checkHighOrderTriangles(GModel *m, std::vector<MElement*> &bad, double &minJGlob)
+static void checkHighOrderTriangles(const char* cc, GModel *m, std::vector<MElement*> &bad, double &minJGlob)
{
bad.clear();
minJGlob = 1.0;
+ double minGGlob = 1.0;
+ double avg = 0.0;
+ int count = 0, nbfair=0;
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); i++){
MTriangle *t = (*it)->triangles[i];
- double disto = t->distoShapeMeasure();
+ double disto_ = t->distoShapeMeasure();
+ double gamma_ = t->gammaShapeMeasure();
+ double disto = disto_;
+ minJGlob = std::min(minJGlob, disto);
+ minGGlob = std::min(minGGlob, gamma_);
+ avg += disto; count++;
+ if (disto < 0) bad.push_back(t);
+ else if (disto < 0.2) nbfair++;
+ }
+ }
+ if (minJGlob > 0) Msg::Info("%s : Worst Triangle Smoothness %g Gamma %g NbFair = %d", cc,minJGlob, minGGlob,nbfair );
+ else Msg::Warning("%s : Worst Triangle Smoothness %g (%d negative jacobians) Worst Gamma %g Avg Smoothness %g",cc, minJGlob, bad.size(),minGGlob,avg/(count?count:1));
+}
+
+
+static void checkHighOrderTetrahedron(const char* cc,GModel *m, std::vector<MElement*> &bad, double &minJGlob)
+{
+ bad.clear();
+ minJGlob = 1.0;
+ double minGGlob = 1.0;
+ double avg = 0.0;
+ int count = 0, nbfair=0;
+ for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
+ for(unsigned int i = 0; i < (*it)->tetrahedra.size(); i++){
+ MTetrahedron *t = (*it)->tetrahedra[i];
+ double disto_ = t->distoShapeMeasure();
+ double gamma_ = t->gammaShapeMeasure();
+ double disto = disto_;
minJGlob = std::min(minJGlob, disto);
+ minGGlob = std::min(minGGlob, gamma_);
+ avg += disto; count++;
if (disto < 0) bad.push_back(t);
+ else if (disto < 0.2) nbfair++;
}
}
- if (minJGlob > 0) Msg::Info("Worst Element Smoothness %g", minJGlob);
- else Msg::Warning("Worst Element Smoothness %g", minJGlob);
+ if (minJGlob > 0) Msg::Info("%s : Worst Tetrahedron Smoothness %g Gamma %g NbFair = %d", cc, minJGlob, minGGlob,nbfair );
+ else Msg::Warning("%s : Worst Tetrahedron Smoothness %g (%d negative jacobians) Worst Gamma %g Avg Smoothness %g", cc, minJGlob, bad.size(),minGGlob,avg/(count?count:1));
}
+
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;
+ const int n = 15;
double D[n][n], X[n][n], Y[n][n], Z[n][n];
FILE *f = fopen(nm,"w");
@@ -1050,14 +1108,11 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
// we do that model face by model face
std::vector<MElement*> bad;
double worst;
+ checkHighOrderTriangles("Before optimization ", m, bad, worst);
if (displ2D){
- checkHighOrderTriangles(m, bad, worst);
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
- if ((*it)->geomType() == GEntity::Plane) {
- gmshSwapHighOrderTriangles(*it);
- if (displ2D) displ2D->smooth(*it);
- }
- // will have to smooth in the planar coordinates, using the metric
+ displ2D->optimize(*it,edgeVertices,faceVertices);
+ checkHighOrderTriangles("After optimization ", m, bad, worst);
}
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
@@ -1065,9 +1120,11 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
displ2D, displ3D);
// smooth the 3D regions
+ checkHighOrderTetrahedron("Before optimization", m, bad, worst);
if (displ3D){
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
displ3D->smooth(*it);
+ checkHighOrderTetrahedron("After optimization", m, bad, worst);
}
if(displ2D){
@@ -1075,8 +1132,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
delete displ3D;
}
- printJacobians(m, "detjIni.pos");
- checkHighOrderTriangles(m, bad, worst);
+ printJacobians(m, "smoothness.pos");
double t2 = Cpu();
Msg::Info("Meshing order %d complete (%g s)", order, t2 - t1);
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index fbdaf071e9..ca6e7d2a3a 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -19,7 +19,20 @@ typedef std::map<std::pair<MVertex*, MVertex*>, std::vector<MVertex*> > edgeCont
// are the high order representation of the face
typedef std::map<MFace, std::vector<MVertex*>, Less_Face> faceContainer;
+// the high order smoother class
+class gmshHighOrderSmoother;
+
void SetOrder1(GModel *m);
void SetOrderN(GModel *m, int order, bool linear=true, bool incomplete=false);
+
+MTriangle* setHighOrder(MTriangle *t,
+ GFace *gf,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ bool linear,
+ bool incomplete,
+ int nPts = 1,
+ gmshHighOrderSmoother *displ2D = 0,
+ gmshHighOrderSmoother *displ3D = 0);
#endif
diff --git a/Mesh/Makefile b/Mesh/Makefile
index fb4a23be92..6c56f66267 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -100,23 +100,24 @@ Generator${OBJEXT}: Generator.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h
Field${OBJEXT}: Field.cpp ../Common/GmshConfig.h \
../contrib/ANN/include/ANN/ANN.h ../Common/Context.h \
../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h Field.h \
- ../Post/PView.h ../Geo/SPoint3.h ../Geo/GeoInterpolation.h ../Geo/Geo.h \
- ../Common/GmshDefines.h ../Geo/gmshSurface.h ../Geo/Pair.h \
+ ../Geo/STensor3.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+ ../Numeric/GmshMatrix.h ../Common/GmshMessage.h ../Numeric/Numeric.h \
+ ../Numeric/GmshMatrix.h ../Post/PView.h ../Geo/GeoInterpolation.h \
+ ../Geo/Geo.h ../Common/GmshDefines.h ../Geo/gmshSurface.h ../Geo/Pair.h \
../Geo/Range.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
- ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
- ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/GmshMessage.h \
- ../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
- ../Common/ListUtils.h ../Geo/SPoint2.h ../Geo/ExtrudeParams.h \
- ../Common/SmoothData.h ../Geo/GModel.h ../Geo/GVertex.h \
- ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
- ../Geo/SBoundingBox3d.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 ../Post/OctreePost.h ../Common/Octree.h \
- ../Common/OctreeInternals.h ../Post/PViewDataList.h ../Post/PViewData.h \
- ../Geo/MVertex.h ../Geo/SPoint2.h ../Geo/SPoint3.h
+ ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Common/ListUtils.h \
+ ../Common/TreeUtils.h ../Common/avl.h ../Common/ListUtils.h \
+ ../Geo/SPoint2.h ../Geo/ExtrudeParams.h ../Common/SmoothData.h \
+ ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
+ ../Geo/SPoint3.h ../Geo/SBoundingBox3d.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 \
+ ../Post/OctreePost.h ../Common/Octree.h ../Common/OctreeInternals.h \
+ ../Post/PViewDataList.h ../Post/PViewData.h ../Geo/MVertex.h \
+ ../Geo/SPoint2.h ../Geo/SPoint3.h
gmshSmoothHighOrder${OBJEXT}: gmshSmoothHighOrder.cpp ../Geo/MLine.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
@@ -189,15 +190,16 @@ meshGFace${OBJEXT}: meshGFace.cpp meshGFace.h meshGFaceBDS.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 ../Geo/MLine.h ../Geo/MElement.h \
- ../Geo/MQuadrangle.h ../Geo/MElement.h ../Common/Context.h \
- ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
+ ../Geo/GEdgeCompound.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/GEdge.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 ../Geo/MLine.h \
+ ../Geo/MElement.h ../Geo/MQuadrangle.h ../Geo/MElement.h \
+ ../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h \
../Numeric/Numeric.h ../Numeric/GmshMatrix.h BDS.h qualityMeasures.h \
- Field.h ../Post/PView.h ../Common/OS.h HighOrder.h
+ Field.h ../Geo/STensor3.h ../Geo/SVector3.h ../Post/PView.h \
+ ../Common/OS.h HighOrder.h
meshGFaceTransfinite${OBJEXT}: meshGFaceTransfinite.cpp meshGFace.h \
../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Geo/GPoint.h \
@@ -249,7 +251,8 @@ meshGFaceBDS${OBJEXT}: meshGFaceBDS.cpp ../Common/GmshMessage.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 ../Numeric/Numeric.h ../Numeric/GmshMatrix.h \
- BDS.h qualityMeasures.h Field.h ../Post/PView.h ../Common/OS.h
+ BDS.h qualityMeasures.h Field.h ../Geo/STensor3.h ../Geo/SVector3.h \
+ ../Post/PView.h ../Common/OS.h
meshGFaceDelaunayInsertion${OBJEXT}: meshGFaceDelaunayInsertion.cpp \
../Common/GmshMessage.h ../Numeric/GmshPredicates.h BackgroundMesh.h \
meshGFaceDelaunayInsertion.h ../Geo/MTriangle.h ../Geo/MElement.h \
@@ -264,7 +267,7 @@ meshGFaceDelaunayInsertion${OBJEXT}: meshGFaceDelaunayInsertion.cpp \
../Geo/GVertex.h ../Geo/GEntity.h ../Geo/GPoint.h ../Geo/SPoint2.h \
../Geo/SVector3.h ../Geo/SPoint3.h ../Geo/SPoint2.h ../Geo/SPoint2.h \
../Geo/SVector3.h ../Geo/Pair.h ../Numeric/Numeric.h \
- ../Numeric/GmshMatrix.h
+ ../Numeric/GmshMatrix.h ../Geo/STensor3.h ../Geo/SVector3.h
meshGFaceOptimize${OBJEXT}: meshGFaceOptimize.cpp meshGFaceOptimize.h \
../Geo/MElement.h ../Common/GmshDefines.h ../Geo/MVertex.h \
../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
@@ -417,7 +420,8 @@ BackgroundMesh${OBJEXT}: BackgroundMesh.cpp ../Common/GmshMessage.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 Field.h ../Post/PView.h
+ ../Geo/SPoint3.h ../Geo/SBoundingBox3d.h Field.h ../Geo/STensor3.h \
+ ../Geo/SVector3.h ../Post/PView.h
qualityMeasures${OBJEXT}: qualityMeasures.cpp qualityMeasures.h BDS.h \
../Common/GmshMessage.h ../Geo/MVertex.h ../Geo/SPoint2.h \
../Geo/SPoint3.h ../Geo/MTriangle.h ../Geo/MElement.h \
@@ -467,16 +471,16 @@ HighOrder${OBJEXT}: HighOrder.cpp HighOrder.h ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/SBoundingBox3d.h ../Geo/MFace.h ../Geo/MVertex.h \
../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/MEdge.h ../Geo/MVertex.h \
../Geo/SVector3.h ../Geo/SVector3.h gmshSmoothHighOrder.h \
+ ../Numeric/GmshMatrix.h ../Common/GmshConfig.h ../Common/GmshMessage.h \
../Geo/MLine.h ../Geo/MElement.h ../Common/GmshDefines.h \
- ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MFace.h ../Common/GmshMessage.h \
- ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h \
- ../Common/GmshConfig.h ../Numeric/Gauss.h ../Geo/MTriangle.h \
- ../Geo/MElement.h ../Geo/MQuadrangle.h ../Geo/MElement.h \
- ../Geo/MTetrahedron.h ../Geo/MElement.h ../Geo/MHexahedron.h \
- ../Geo/MElement.h ../Geo/MPrism.h ../Geo/MElement.h ../Geo/MPyramid.h \
- ../Geo/MElement.h ../Common/OS.h ../Numeric/Numeric.h \
- ../Numeric/GmshMatrix.h ../Common/Context.h ../Geo/CGNSOptions.h \
- ../Mesh/meshPartitionOptions.h
+ ../Geo/MVertex.h ../Geo/MEdge.h ../Geo/MFace.h \
+ ../Numeric/FunctionSpace.h ../Numeric/GmshMatrix.h ../Numeric/Gauss.h \
+ ../Geo/MTriangle.h ../Geo/MElement.h ../Geo/MQuadrangle.h \
+ ../Geo/MElement.h ../Geo/MTetrahedron.h ../Geo/MElement.h \
+ ../Geo/MHexahedron.h ../Geo/MElement.h ../Geo/MPrism.h \
+ ../Geo/MElement.h ../Geo/MPyramid.h ../Geo/MElement.h ../Common/OS.h \
+ ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Common/Context.h \
+ ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
meshPartition${OBJEXT}: meshPartition.cpp ../Common/GmshConfig.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 \
diff --git a/Mesh/gmshSmoothHighOrder.cpp b/Mesh/gmshSmoothHighOrder.cpp
index 1967a77eac..4ec085e0db 100644
--- a/Mesh/gmshSmoothHighOrder.cpp
+++ b/Mesh/gmshSmoothHighOrder.cpp
@@ -28,6 +28,248 @@
#define SQU(a) ((a)*(a))
+int optimalLocationPN_ (GFace *gf, const MEdge &me, MTriangle *t1, MTriangle *t2,gmshHighOrderSmoother *s);
+static int _gmshSwapHighOrderTriangles(GFace *gf);
+static int _gmshSwapHighOrderTriangles(GFace *gf,edgeContainer&,faceContainer&, gmshHighOrderSmoother *s);
+static int _gmshFindOptimalLocationsP2(GFace *gf, gmshHighOrderSmoother *s);
+static int _gmshFindOptimalLocationsPN(GFace *gf,gmshHighOrderSmoother *s);
+
+static double shapeMeasure (MElement *e) {
+ const double d1 = e->distoShapeMeasure();
+ const double d2 = e->gammaShapeMeasure();
+ return d1;
+}
+
+void gmshHighOrderSmoother::moveTo (MVertex *v, const std::map<MVertex*,SVector3> &m) const
+{
+ std::map<MVertex*,SVector3>::const_iterator it = m.find(v);
+ if (it != m.end()){
+ v->x() = it->second.x();
+ v->y() = it->second.y();
+ v->z() = it->second.z();
+ }
+}
+
+void gmshHighOrderSmoother::moveToStraightSidedLocation (MVertex *v) const
+{
+ moveTo ( v, _straightSidedLocation);
+}
+
+void gmshHighOrderSmoother::moveToTargetLocation (MVertex *v) const
+{
+ moveTo ( v, _targetLocation);
+}
+
+void gmshHighOrderSmoother::moveToStraightSidedLocation (MElement *e) const
+{
+ for (int i=0;i<e->getNumVertices();i++)
+ moveToStraightSidedLocation(e->getVertex(i));
+}
+
+void gmshHighOrderSmoother::moveToTargetLocation (MElement *e) const
+{
+ for (int i=0;i<e->getNumVertices();i++)
+ moveToTargetLocation(e->getVertex(i));
+}
+
+void gmshHighOrderSmoother::updateTargetLocation (MVertex*v, const SPoint3 &p3, const SPoint2 &p2)
+{
+ v->x() = p3.x();
+ v->y() = p3.y();
+ v->z() = p3.z();
+ v->setParameter(0,p2.x());
+ v->setParameter(1,p2.y());
+ _targetLocation[v] = SVector3(p3.x(),p3.y(),p3.z());
+}
+
+
+struct p2data{
+ GFace *gf;
+ MTriangle *t1,*t2;
+ MVertex *n12;
+ gmshMatrix<double> *m1,*m2;
+ gmshHighOrderSmoother *s;
+ p2data (GFace *_gf,MTriangle *_t1,MTriangle *_t2,MVertex *_n12,gmshHighOrderSmoother *_s)
+ : gf(_gf),t1(_t1),t2(_t2),n12(_n12),s(_s){
+ gmshElasticityTerm el(0,1.e3,.3333);
+ s->moveToStraightSidedLocation(t1);
+ s->moveToStraightSidedLocation(t2);
+ m1 = new gmshMatrix<double>(3*t1->getNumVertices(),
+ 3*t1->getNumVertices());
+ m2 = new gmshMatrix<double>(3*t2->getNumVertices(),
+ 3*t2->getNumVertices());
+ el.elementMatrix(t1,*m1);
+ el.elementMatrix(t2,*m2);
+ s->moveToTargetLocation(t1);
+ s->moveToTargetLocation(t2);
+ }
+ ~p2data(){
+ delete m1;
+ delete m2;
+ }
+};
+
+struct pNdata{
+ GFace *gf;
+ MTriangle *t1,*t2;
+ const std::vector<MVertex*> &n;
+ gmshMatrix<double> *m1,*m2;
+ gmshHighOrderSmoother *s;
+ pNdata (GFace *_gf,MTriangle *_t1,MTriangle *_t2,const std::vector<MVertex*> &_n,gmshHighOrderSmoother *_s)
+ : gf(_gf),t1(_t1),t2(_t2),n(_n),s(_s)
+ {
+ gmshElasticityTerm el(0,1.e3,.3333);
+ s->moveToStraightSidedLocation(t1);
+ s->moveToStraightSidedLocation(t2);
+ m1 = new gmshMatrix<double>(3*t1->getNumVertices(),
+ 3*t1->getNumVertices());
+ m2 = new gmshMatrix<double>(3*t2->getNumVertices(),
+ 3*t2->getNumVertices());
+ el.elementMatrix(t1,*m1);
+ el.elementMatrix(t2,*m2);
+ s->moveToTargetLocation(t1);
+ s->moveToTargetLocation(t2);
+ }
+ ~pNdata(){
+ delete m1;
+ delete m2;
+ }
+};
+
+
+static double _DeformationEnergy (MElement *e,
+ gmshMatrix<double> *K,
+ gmshHighOrderSmoother *s){
+
+ int N = e->getNumVertices();
+ gmshVector<double> Kdx(N*3),dx(N*3);
+
+ dx.scale(0.0);
+ Kdx.scale(0.0);
+ int C = 0;
+ for (int i=0;i<N;i++){
+ SVector3 disp = s->getDisplacement(e->getVertex(i));
+ SVector3 str = s->getSSL(e->getVertex(i));
+ dx(i) = disp.x();
+ dx(i+N) = disp.y();
+ dx(i+2*N) = disp.z();
+ if (0 && (fabs(disp.x())>1.e-12 ||fabs(disp.y())>1.e-12))
+ printf("%6d (%12.5E %12.5E %12.5E -- %12.5E %12.5E %12.5E -- %12.5E %12.5E %12.5E)\n",
+ e->getVertex(i)->getNum(),
+ disp.x(),disp.y(),disp.z(),
+ str.x(),str.y(),str.z(),e->getVertex(i)->x(),e->getVertex(i)->y(),
+ e->getVertex(i)->z());
+ }
+
+ K->mult(dx,Kdx);
+ const double E = Kdx * dx;
+ return E;
+}
+
+static double _DeformationEnergy (p2data *p){
+
+ // printf("DEFEN \n");
+
+ return
+ _DeformationEnergy (p->t1,p->m1,p->s) +
+ _DeformationEnergy (p->t2,p->m2,p->s);
+}
+
+static double _DeformationEnergy (pNdata *p){
+
+ return _DeformationEnergy (p->t1,p->m1,p->s) +
+ _DeformationEnergy (p->t2,p->m2,p->s);
+}
+
+
+static double _function_p2tB (gmshVector<double> &x, void *data){
+ p2data *p = (p2data*)data;
+ GPoint gp12 = p->gf->point(SPoint2(x(0),x(1)));
+ double xx = p->n12->x();
+ double yy = p->n12->y();
+ double zz = p->n12->z();
+ p->n12->x() = gp12.x();
+ p->n12->y() = gp12.y();
+ p->n12->z() = gp12.z();
+
+ double E = _DeformationEnergy(p);
+ // printf("%12.5E %12.5E -- %12.5E %12.5E => %12.5E %12.5E E %12.5E\n",x(0),x(1),xx,yy,gp12.x(),gp12.y(),E);
+
+ p->n12->x() = xx;
+ p->n12->y() = yy;
+ p->n12->z() = zz;
+ return E;
+}
+
+
+
+static double _function_p2t (gmshVector<double> &x, void *data){
+ p2data *p = (p2data*)data;
+ GPoint gp12 = p->gf->point(SPoint2(x(0),x(1)));
+ // printf("%g %g = %g %g\n",x(0),x(1),gp12.x(),gp12.y());
+ double xx = p->n12->x();
+ double yy = p->n12->y();
+ double zz = p->n12->z();
+ p->n12->x() = gp12.x();
+ p->n12->y() = gp12.y();
+ p->n12->z() = gp12.z();
+ double q = std::min(shapeMeasure(p->t1),shapeMeasure(p->t2));
+ p->n12->x() = xx;
+ p->n12->y() = yy;
+ p->n12->z() = zz;
+ return -q;
+}
+
+static double _function_pNt (gmshVector<double> &x, void *data){
+ pNdata *p = (pNdata*)data;
+ static double xx[256];
+ static double yy[256];
+ static double zz[256];
+ for (int i=0;i<p->n.size();i++){
+ GPoint gp12 = p->gf->point(SPoint2(x(2*i),x(2*i+1)));
+ // printf("%g %g = %g %g\n",x(0),x(1),gp12.x(),gp12.y());
+ xx[i] = p->n[i]->x();
+ yy[i] = p->n[i]->y();
+ zz[i] = p->n[i]->z();
+ p->n[i]->x() = gp12.x();
+ p->n[i]->y() = gp12.y();
+ p->n[i]->z() = gp12.z();
+ }
+ double q = std::min(shapeMeasure(p->t1),shapeMeasure(p->t2));
+ for (int i=0;i<p->n.size();i++){
+ p->n[i]->x() = xx[i];
+ p->n[i]->y() = yy[i];
+ p->n[i]->z() = zz[i];
+ }
+ return -q;
+}
+
+static double _function_pNtB (gmshVector<double> &x, void *data){
+ pNdata *p = (pNdata*)data;
+ static double xx[256];
+ static double yy[256];
+ static double zz[256];
+ for (int i=0;i<p->n.size();i++){
+ GPoint gp12 = p->gf->point(SPoint2(x(2*i),x(2*i+1)));
+ // printf("%g %g = %g %g\n",x(0),x(1),gp12.x(),gp12.y());
+ xx[i] = p->n[i]->x();
+ yy[i] = p->n[i]->y();
+ zz[i] = p->n[i]->z();
+ p->n[i]->x() = gp12.x();
+ p->n[i]->y() = gp12.y();
+ p->n[i]->z() = gp12.z();
+ }
+ double E = _DeformationEnergy(p);
+ for (int i=0;i<p->n.size();i++){
+ p->n[i]->x() = xx[i];
+ p->n[i]->y() = yy[i];
+ p->n[i]->z() = zz[i];
+ }
+ // printf("E(%g,%g) = %g\n",x(0),x(1),E);
+ return E;
+}
+
+
void getDistordedElements(const std::vector<MElement*> & v,
const double & threshold,
std::vector<MElement*> & d,
@@ -77,7 +319,7 @@ void addOneLayer(const std::vector<MElement*> & v,
}
}
-void gmshHighOrderSmoother::smooth(GFace *gf)
+void gmshHighOrderSmoother::smooth(GFace *gf, bool metric)
{
std::vector<MElement*> v;
@@ -85,7 +327,8 @@ void gmshHighOrderSmoother::smooth(GFace *gf)
v.insert(v.begin(), gf->quadrangles.begin(),gf->quadrangles.end());
Msg::Info("Smoothing high order mesh : model face %d (%d elements)", gf->tag(),
v.size());
- smooth(v);
+ if (metric)smooth_metric(v,gf);
+ else smooth(v);
}
void gmshHighOrderSmoother::smooth(GRegion *gr)
@@ -99,19 +342,106 @@ void gmshHighOrderSmoother::smooth(GRegion *gr)
smooth(v);
}
-void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
+// Use elastic solver to move the nodes
+
+// compute the stiffness matrix of an element
+// that correspond to the deformation of a
+// straight sided element to a curvilinear one
+
+
+
+
+
+void gmshHighOrderSmoother::optimize(GFace * gf,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices){
+
+ // if (gf->geomType() != GEntity::Plane)return;
+
+ while (1) {
+ // relocate the vertices using elliptic smoother
+ // smooth(gf);
+ // for (int i=0;i<20;i++){
+ // _gmshFindOptimalLocationsPN(gf,this);
+ // _gmshFindOptimalLocationsPN(gf,this);
+ // }
+ // then try to swap for better configurations
+
+ smooth(gf,true);
+ int nbSwap = _gmshSwapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
+ //smooth(gf,true);
+ // smooth(gf,true);
+ // smooth(gf,true);
+ // smooth(gf,true);
+ break;
+
+ // printf("%d %g\n",nbSwap,_MIDDLE);
+ // if the smoothing procedure has been able to
+ // move the nodes to their right location, break !
+ // break;
+ if (_MIDDLE == 1.0) break;
+ // if (!nbSwap){
+ // printf("Cannot move thet nodes to their optimal location, splits have to be added\n");
+ // break;
+ // }
+ }
+}
+
+
+void gmshHighOrderSmoother::computeMetricVector (GFace *gf,
+ MElement *e,
+ gmshMatrix<double> &J,
+ gmshMatrix<double> &JT,
+ gmshVector<double> &D){
+ int nbNodes = e->getNumVertices();
+ for (int j = 0; j < nbNodes; j++){
+ SPoint2 param;
+ reparamMeshVertexOnFace(e->getVertex(j), gf, param);
+ Pair<SVector3,SVector3> der = gf->firstDer(param);
+ int XJ = j;
+ int YJ = j+nbNodes;
+ int ZJ = j+2*nbNodes;
+ int UJ = j;
+ int VJ = j+nbNodes;
+ J(XJ,UJ) = der.first().x();
+ J(YJ,UJ) = der.first().y();
+ J(ZJ,UJ) = der.first().z();
+ J(XJ,VJ) = der.second().x();
+ J(YJ,VJ) = der.second().y();
+ J(ZJ,VJ) = der.second().z();
+
+ JT(UJ,XJ) = der.first().x();
+ JT(UJ,YJ) = der.first().y();
+ JT(UJ,ZJ) = der.first().z();
+ JT(VJ,XJ) = der.second().x();
+ JT(VJ,YJ) = der.second().y();
+ JT(VJ,ZJ) = der.second().z();
+
+ GPoint gp = gf->point(param);
+ SVector3 ss = getSSL(e->getVertex(j));
+ D(XJ) = gp.x() - ss.x();
+ D(YJ) = gp.y() - ss.y();
+ D(ZJ) = gp.z() - ss.z();
+ }
+}
+
+
+void gmshHighOrderSmoother::smooth_metric ( std::vector<MElement*> & all, GFace *gf)
{
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
+ // lsys->setNoisy(1);
+ lsys->setGmres(1);
+ lsys->setPrec(1.e-7);
gmshAssembler<double> myAssembler(lsys);
gmshElasticityTerm El(0, 1.0, .333, getTag());
- std::vector<MElement*> v, layer;
+ std::vector<MElement*> layer,v;
double minD;
- getDistordedElements ( all, 0.5, v,minD);
+ getDistordedElements ( all, .6, v,minD);
- // printf("%d elements / %d distorted min Disto = %g\n",all.size(),v.size(), minD);
+ Msg::Debug("%d elements / %d distorted min Disto = %g",all.size(),v.size(), minD);
if (!v.size()) return;
@@ -122,24 +452,212 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
}
// 3 -> .4
- printf("%d elements after adding %d layers\n", (int)v.size(), nbLayers);
+ Msg::Debug("%d elements after adding %d layers", (int)v.size(), nbLayers);
addOneLayer(all, v, layer);
// printf("%d elements in the next layer\n",layer.size());
+ std::set<MVertex*>::iterator it;
+ std::map<MVertex*,SVector3>::iterator its;
+ std::map<MVertex*,SVector3>::iterator itpresent;
+ std::map<MVertex*,SVector3>::iterator ittarget;
+ std::set<MVertex*> verticesToMove;
+
for (unsigned int i = 0; i < layer.size(); i++){
for (int j = 0; j < layer[i]->getNumVertices(); j++){
MVertex *vert = layer[i]->getVertex(j);
myAssembler.fixVertex(vert, 0, getTag(), 0);
myAssembler.fixVertex(vert, 1, getTag(), 0);
- myAssembler.fixVertex(vert, 2, getTag(), 0);
}
}
+
+ // printf("%d vertices \n", _displ.size());
+
+ for (unsigned int i = 0; i < v.size(); i++){
+ for (int j = 0; j < v[i]->getNumVertices(); j++){
+ MVertex *vert = v[i]->getVertex(j);
+ // printf("%d %d %d v\n",i,j,v[i]->getNumVertices());
+ its = _straightSidedLocation.find(vert);
+ if (its == _straightSidedLocation.end()){
+ _straightSidedLocation[vert] =
+ SVector3(vert->x(),vert->y(),vert->z());
+ _targetLocation[vert] =
+ SVector3(vert->x(),vert->y(),vert->z());
+ }
+ else{
+ vert->x() = its->second.x();
+ vert->y() = its->second.y();
+ vert->z() = its->second.z();
+ if (vert->onWhat()->dim() < _dim){
+ myAssembler.fixVertex ( vert , 0 , getTag() , 0);
+ myAssembler.fixVertex ( vert , 1 , getTag() , 0);
+ }
+ }
+ }
+ }
+
+ // number the other DOFs
+ for (unsigned int i = 0; i < v.size(); i++){
+ for (int j = 0; j < v[i]->getNumVertices(); j++){
+ MVertex *vert = v[i]->getVertex(j);
+ myAssembler.numberVertex(vert, 0, getTag());
+ myAssembler.numberVertex(vert, 1, getTag());
+ verticesToMove.insert(vert);
+ }
+ }
+
+
+ double dx0 = smooth_metric_ ( v, gf, myAssembler, verticesToMove,El);
+ double dx = dx0;
+ // printf(" dx0 = %12.5E\n",dx0);
+ int iter = 0;
+ while(1){
+ double dx2 = smooth_metric_ ( v,gf, myAssembler, verticesToMove,El);
+ // printf(" dx2 = %12.5E\n",dx2);
+ if (fabs(dx2-dx) < 1.e-4 * dx0)break;
+ if (iter++ > 10)break;
+ dx = dx2;
+ }
+
+ for (it = verticesToMove.begin() ; it != verticesToMove.end() ; ++it){
+ SPoint2 param;
+ reparamMeshVertexOnFace(*it, gf, param);
+ GPoint gp = gf->point(param);
+ if ((*it)->onWhat()->dim() == 2){
+ (*it)->x() = gp.x();
+ (*it)->y() = gp.y();
+ (*it)->z() = gp.z();
+ _targetLocation[*it] = SVector3(gp.x(),gp.y(),gp.z());
+ }
+ else{
+ SVector3 p = _targetLocation[(*it)];
+ (*it)->x() = p.x();
+ (*it)->y() = p.y();
+ (*it)->z() = p.z();
+ }
+ }
+ delete lsys;
+}
+
+
+double gmshHighOrderSmoother::smooth_metric_ ( std::vector<MElement*> & v,
+ GFace *gf,
+ gmshAssembler<double> &myAssembler,
+ std::set<MVertex*> &verticesToMove,
+ gmshElasticityTerm &El)
+{
+ // Msg::Info("%d vertices FIXED %d NUMBERED", myAssembler.sizeOfF()
+ // , myAssembler.sizeOfR());
+
+ std::set<MVertex*>::iterator it;
+
+ if (myAssembler.sizeOfR()){
+
+ for (int i=0;i<v.size();i++){
+ MElement *e = v[i];
+ int nbNodes = e->getNumVertices();
+ const int n2 = 2*nbNodes;
+ const int n3 = 3*nbNodes;
+
+ gmshMatrix<double> K33 (n3,n3);
+ gmshMatrix<double> K22 (n2,n2);
+ gmshMatrix<double> J32 (n3,n2);
+ gmshMatrix<double> J23 (n2,n3);
+ gmshVector<double> D3 (n3);
+ gmshVector<double> R2 (n2);
+ gmshMatrix<double> J23K33 (n2,n3);
+ K33.set_all(0.0);
+ El.elementMatrix(e,K33);
+ computeMetricVector (gf,e,J32,J23,D3);
+ J23K33.gemm(J23,K33,1,0);
+ K22.gemm(J23K33,J32,1,0);
+ J23K33.mult(D3,R2);
+ for (int j = 0; j < n2; j++){
+ MVertex *vR;
+ int iCompR, iFieldR;
+ El.getLocalDofR(e, j, &vR, &iCompR, &iFieldR);
+ myAssembler.assemble(vR, iCompR, iFieldR,-R2(j));
+ for (int k = 0; k < n2; k++){
+ MVertex *vC;
+ int iCompC, iFieldC;
+ El.getLocalDofC(e, k, &vC, &iCompC, &iFieldC);
+ myAssembler.assemble(vR, iCompR, iFieldR,
+ vC, iCompC, iFieldC,
+ K22(j, k));
+ }
+ }
+ }
+ // solve the system
+ myAssembler.systemSolve();
+ }
+
+
+ double dx = 0.0;
+ for (it = verticesToMove.begin() ; it != verticesToMove.end() ; ++it){
+ if ((*it)->onWhat()->dim() == 2){
+ SPoint2 param;
+ reparamMeshVertexOnFace((*it), gf, param);
+ SPoint2 dparam (myAssembler.getDofValue ((*it), 0 ,getTag()),
+ myAssembler.getDofValue ((*it), 1 ,getTag()));
+ SPoint2 newp = param+dparam;
+ dx += newp.x()*newp.x()+newp.y()*newp.y();
+ (*it)->setParameter(0,newp.x());
+ (*it)->setParameter(1,newp.y());
+ }
+ }
+
+ myAssembler.systemClear();
+
+ return dx;
+}
+
+
+void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
+{
+ gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
+ lsys->setNoisy(1);
+ lsys->setGmres(1);
+ lsys->setPrec(1.e-5);
+ gmshAssembler<double> myAssembler(lsys);
+ gmshElasticityTerm El(0, 1.0, .333, getTag());
+ std::vector<MElement*> layer,v;
+
+ double minD;
+
+ getDistordedElements ( all, .5, v,minD);
+
+ Msg::Debug("%d elements / %d distorted min Disto = %g\n",all.size(),v.size(), minD);
+
+ if (!v.size()) return;
+
+ const int nbLayers = 2;
+ for (int i=0;i<nbLayers;i++){
+ addOneLayer ( all, v, layer);
+ v.insert(v.end(),layer.begin(),layer.end());
+ }
+
+ // 3 -> .4
+ Msg::Debug("%d elements after adding %d layers\n", (int)v.size(), nbLayers);
+
+ addOneLayer(all, v, layer);
+
+ // printf("%d elements in the next layer\n",layer.size());
+
+ for (unsigned int i = 0; i < layer.size(); i++){
+ for (int j = 0; j < layer[i]->getNumVertices(); j++){
+ MVertex *vert = layer[i]->getVertex(j);
+ myAssembler.fixVertex(vert, 0, getTag(), 0);
+ myAssembler.fixVertex(vert, 1, getTag(), 0);
+ myAssembler.fixVertex(vert, 2, getTag(), 0);
+ }
+ }
std::map<MVertex*,SVector3>::iterator it;
+ std::map<MVertex*,SVector3>::iterator its;
+ std::map<MVertex*,SVector3>::iterator itt;
std::map<MVertex*,SVector3> verticesToMove;
// printf("%d vertices \n", _displ.size());
@@ -148,13 +666,12 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
for (int j = 0; j < v[i]->getNumVertices(); j++){
MVertex *vert = v[i]->getVertex(j);
// printf("%d %d %d v\n",i,j,v[i]->getNumVertices());
- it = _displ.find(vert);
- if (it != _displ.end()){
- myAssembler.fixVertex ( vert , 0 , getTag() , vert->x()-it->second.x());
- myAssembler.fixVertex ( vert , 1 , getTag() , vert->y()-it->second.y());
- myAssembler.fixVertex ( vert , 2 , getTag() , vert->z()-it->second.z());
- // printf("%g %g vs %g %g\n", it->second.x(), it->second.y(),vert->x(),vert->y());
- verticesToMove[vert] = it->second;
+ its = _straightSidedLocation.find(vert);
+ itt = _targetLocation.find(vert);
+ if (its != _straightSidedLocation.end() && vert->onWhat()->dim() < _dim){
+ myAssembler.fixVertex ( vert , 0 , getTag() , vert->x()-its->second.x());
+ myAssembler.fixVertex ( vert , 1 , getTag() , vert->y()-its->second.y());
+ myAssembler.fixVertex ( vert , 2 , getTag() , vert->z()-its->second.z());
}
// ensure we do not touch any vertex that is on the boundary
else if (vert->onWhat()->dim() < _dim){
@@ -164,14 +681,9 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
}
}
}
-
- // move back high order nodes to their straight sided
- // location
- for (it = verticesToMove.begin() ; it != verticesToMove.end() ; ++it){
- it->first->x() = it->second.x();
- it->first->y() = it->second.y();
- it->first->z() = it->second.z();
- }
+
+ for (unsigned int i = 0; i < v.size(); i++)
+ moveToStraightSidedLocation (v[i]);
// number the other DOFs
for (unsigned int i = 0; i < v.size(); i++){
@@ -197,24 +709,15 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
// solve the system
lsys->systemSolve();
}
-
- // move the nodes that were involved in the process
- // to their new lcation
+
for (it = verticesToMove.begin() ; it != verticesToMove.end() ; ++it){
- it->first->x() = it->first->x() + myAssembler.getDofValue (it->first, 0 ,getTag());
- it->first->y() = it->first->y() + myAssembler.getDofValue (it->first, 1 ,getTag());
- it->first->z() = it->first->z() + myAssembler.getDofValue (it->first, 2 ,getTag());
- // printf("%g %g\n", myAssembler.getDofValue (it->first, 0 ,getTag()), myAssembler.getDofValue (it->first, 1 ,getTag()));
+ it->first->x() += myAssembler.getDofValue (it->first, 0 ,getTag());
+ it->first->y() += myAssembler.getDofValue (it->first, 1 ,getTag());
+ it->first->z() += myAssembler.getDofValue (it->first, 2 ,getTag());
}
// delete matrices and vectors
-
- double minD2;
- getDistordedElements ( v, 0.5, layer,minD2);
-
- printf("Smooting efficiency %g -> %g\n",minD,minD2);
-
-
+
delete lsys;
}
@@ -232,6 +735,120 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
n1 n14 n4
*/
+/*
+ n34 is the new vertex, we'd like to put it at the
+ best place i.e. at a place where it maximizes the
+ minimal smoothness of the neighboring elements.
+
+ One can actually maximize this value or one can use
+ a standard interpolation scheme (transfinite or elliptic)
+ to place the point.
+*/
+
+
+static void getNodesP2 (const MEdge &me, MTriangle *t1, MTriangle *t2,
+ MVertex * &n1,MVertex * &n2,MVertex * &n3,MVertex * &n4,
+ MVertex * &n12,MVertex * &n14,MVertex * &n24,MVertex * &n23,MVertex * &n13){
+
+ n1 = me.getVertex(0);
+ n2 = me.getVertex(1);
+
+ if (t1->getVertex(0) != n1 && t1->getVertex(0) != n2) n3 = t1->getVertex(0);
+ else if (t1->getVertex(1) != n1 && t1->getVertex(1) != n2) n3 = t1->getVertex(1);
+ else if (t1->getVertex(2) != n1 && t1->getVertex(2) != n2) n3 = t1->getVertex(2);
+ int iStart = 0;
+ for (;iStart<3;iStart++)
+ if (t1->getVertex(iStart) == n1)
+ break;
+ if (n2 == t1->getVertex((iStart+1)%3)){
+ n12 = t1->getVertex((iStart+0)%3 + 3);
+ n23 = t1->getVertex((iStart+1)%3 + 3);
+ n13 = t1->getVertex((iStart+2)%3 + 3);
+ }
+ else{
+ n13 = t1->getVertex((iStart+0)%3 + 3);
+ n23 = t1->getVertex((iStart+1)%3 + 3);
+ n12 = t1->getVertex((iStart+2)%3 + 3);
+ }
+
+ if (t2->getVertex(0) != n1 && t2->getVertex(0) != n2) n4 = t2->getVertex(0);
+ else if (t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
+ else if (t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
+ iStart = 0;
+ for (;iStart<3;iStart++)
+ if (t2->getVertex(iStart) == n1)
+ break;
+ if (n2 == t2->getVertex((iStart+1)%3)){
+ n24 = t2->getVertex((iStart+1)%3 + 3);
+ n14 = t2->getVertex((iStart+2)%3 + 3);
+ }
+ else{
+ n14 = t2->getVertex((iStart+0)%3 + 3);
+ n24 = t2->getVertex((iStart+1)%3 + 3);
+ }
+}
+
+static void getNodesPN (const MEdge &me, MTriangle *t1, MTriangle *t2,
+ MVertex * &n1,MVertex * &n2,MVertex * &n3,MVertex * &n4,
+ std::vector<MVertex *> &n12,
+ std::vector<MVertex *> &n14,
+ std::vector<MVertex *> &n24,
+ std::vector<MVertex *> &n23,
+ std::vector<MVertex *> &n13){
+
+ n1 = me.getVertex(0);
+ n2 = me.getVertex(1);
+
+ int order = t1->getPolynomialOrder();
+
+ if (t1->getVertex(0) != n1 && t1->getVertex(0) != n2) n3 = t1->getVertex(0);
+ else if (t1->getVertex(1) != n1 && t1->getVertex(1) != n2) n3 = t1->getVertex(1);
+ else if (t1->getVertex(2) != n1 && t1->getVertex(2) != n2) n3 = t1->getVertex(2);
+ int iStart = 0;
+ for (;iStart<3;iStart++)
+ if (t1->getVertex(iStart) == n1)
+ break;
+ if (n2 == t1->getVertex((iStart+1)%3)){
+ int start12 = 3 + ((iStart+0)%3) * (order-1);
+ for (int i=0;i<order-1;i++)n12.push_back(t1->getVertex(start12+i));
+ int start23 = 3 + ((iStart+1)%3) * (order-1);
+ for (int i=0;i<order-1;i++)n23.push_back(t1->getVertex(start23+i));
+ int start13 = 3 + ((iStart+2)%3) * (order-1);
+ for (int i=0;i<order-1;i++)n13.push_back(t1->getVertex(start13+i));
+ }
+ else{
+ int start13 = 3 + ((iStart+0)%3) * (order-1);
+ for (int i=order-2;i>=0;i--)n13.push_back(t1->getVertex(start13+i));
+ int start23 = 3 + ((iStart+1)%3) * (order-1);
+ for (int i=order-2;i>=0;i--)n23.push_back(t1->getVertex(start23+i));
+ int start12 = 3 + ((iStart+2)%3) * (order-1);
+ for (int i=order-2;i>=0;i--)n12.push_back(t1->getVertex(start12+i));
+ }
+
+ if (t2->getVertex(0) != n1 && t2->getVertex(0) != n2) n4 = t2->getVertex(0);
+ else if (t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
+ else if (t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
+ iStart = 0;
+ for (;iStart<3;iStart++)
+ if (t2->getVertex(iStart) == n1)
+ break;
+
+ if (n2 == t2->getVertex((iStart+1)%3)){
+ int start24 = 3 + ((iStart+1)%3) * (order-1);
+ for (int i=0;i<order-1;i++)n24.push_back(t2->getVertex(start24+i));
+ int start14 = 3 + ((iStart+2)%3) * (order-1);
+ for (int i=0;i<order-1;i++)n14.push_back(t2->getVertex(start14+i));
+ }
+ else{
+ int start14 = 3 + ((iStart+0)%3) * (order-1);
+ for (int i=order-2;i>=0;i--)n14.push_back(t1->getVertex(start14+i));
+ int start24 = 3 + ((iStart+1)%3) * (order-1);
+ for (int i=order-2;i>=0;i--)n24.push_back(t1->getVertex(start24+i));
+ }
+}
+
+
+
static double surfaceTriangleUV(SPoint2 &v1, SPoint2 &v2, SPoint2 &v3)
{
const double v12[2] = {v2.x() - v1.x(),v2.y() - v1.y()};
@@ -239,8 +856,6 @@ static double surfaceTriangleUV(SPoint2 &v1, SPoint2 &v2, SPoint2 &v3)
return 0.5 * fabs (v12[0] * v13[1] - v12[1] * v13[0]);
}
-
-
struct swap_triangles_p2
{
MTriangle *t1, *t2;
@@ -251,61 +866,66 @@ struct swap_triangles_p2
MVertex *n1, *n2, *n3, *n4;
MVertex *n12, *n13, *n23, *n14, *n24;
MVertex *n34;
+
+
+ double evalConfiguration (GFace *gf, SPoint2 & p) const {
+ GPoint gp34 = gf->point(p);
+ MFaceVertex _test (gp34.x(),gp34.y(),gp34.z(),gf,p.x(),p.y());
+ std::vector<MVertex *>vv;
+ vv.push_back(n13);vv.push_back(&_test);vv.push_back(n14);
+ MTriangleN t3_test (n1,n3,n4,vv,2,t1->getNum(),t1->getPartition());
+ vv.clear();
+ vv.push_back(&_test);vv.push_back(n23);vv.push_back(n24);
+ MTriangleN t4_test (n4,n3,n2,vv,2,t2->getNum(),t2->getPartition());
+ const double q3 = shapeMeasure(&t3_test);
+ const double q4 = shapeMeasure(&t4_test);
+ return std::min(q3,q4);
+ }
+
+ MVertex *optimalLocation (GFace *gf,
+ SPoint2 &p3,
+ SPoint2 &p4) const {
+ SPoint2 p34_linear = (p3+p4)*.5;
+
+
+ GPoint gp34 = gf->point(p34_linear);
+ MFaceVertex *_test = new MFaceVertex (gp34.x(),gp34.y(),gp34.z(),
+ gf,p34_linear.x(),p34_linear.y());
+ std::vector<MVertex *>vv;
+ vv.push_back(n13);vv.push_back(_test);vv.push_back(n14);
+ MTriangleN t3_test (n1,n3,n4,vv,2,t1->getNum(),t1->getPartition());
+ vv.clear();
+ vv.push_back(_test);vv.push_back(n23);vv.push_back(n24);
+ MTriangleN t4_test (n4,n3,n2,vv,2,t2->getNum(),t2->getPartition());
+ p2data data(gf,&t3_test,&t4_test,_test,0);
+ gmshVector<double> pp(2);
+ pp(0) = p34_linear.x();
+ pp(1) = p34_linear.y();
+ double opti = minimize_grad_fd (_function_p2tB, pp, &data);
+ return _test;
+ }
+
+
swap_triangles_p2(const MEdge &me, MTriangle *_t1, MTriangle *_t2, GFace *gf)
: t1(_t1), t2(_t2),n1(0),n2(0),n3(0),n12(0),n13(0), n23(0), n14(0), n24(0),n34(0)
{
- n1 = me.getVertex(0);
- n2 = me.getVertex(1);
-
- if (t1->getVertex(0) != n1 && t1->getVertex(0) != n2) n3 = t1->getVertex(0);
- else if (t1->getVertex(1) != n1 && t1->getVertex(1) != n2) n3 = t1->getVertex(1);
- else if (t1->getVertex(2) != n1 && t1->getVertex(2) != n2) n3 = t1->getVertex(2);
- int iStart = 0;
- for (;iStart<3;iStart++)
- if (t1->getVertex(iStart) == n1)
- break;
- if (n2 == t1->getVertex((iStart+1)%3)){
- n12 = t1->getVertex((iStart+0)%3 + 3);
- n23 = t1->getVertex((iStart+1)%3 + 3);
- n13 = t1->getVertex((iStart+2)%3 + 3);
- }
- else{
- n13 = t1->getVertex((iStart+0)%3 + 3);
- n23 = t1->getVertex((iStart+1)%3 + 3);
- n12 = t1->getVertex((iStart+2)%3 + 3);
- }
-
- if (t2->getVertex(0) != n1 && t2->getVertex(0) != n2) n4 = t2->getVertex(0);
- else if (t2->getVertex(1) != n1 && t2->getVertex(1) != n2) n4 = t2->getVertex(1);
- else if (t2->getVertex(2) != n1 && t2->getVertex(2) != n2) n4 = t2->getVertex(2);
- iStart = 0;
- for (;iStart<3;iStart++)
- if (t2->getVertex(iStart) == n1)
- break;
- if (n2 == t2->getVertex((iStart+1)%3)){
- n24 = t2->getVertex((iStart+1)%3 + 3);
- n14 = t2->getVertex((iStart+2)%3 + 3);
- }
- else{
- n14 = t2->getVertex((iStart+0)%3 + 3);
- n24 = t2->getVertex((iStart+1)%3 + 3);
- }
- //printf("%p %p %p %p %p %p %p %p %p\n",n1,n2,n3,n4,n12,n23,n13,n24,n14);
+ const double qold1 = shapeMeasure(t1);
+ const double qold2 = shapeMeasure(t2);
- SPoint2 p1,p2,p3,p4;
+ //printf("%p %p %p %p %p %p %p %p %p\n",n1,n2,n3,n4,n12,n23,n13,n24,n14);
+ getNodesP2 (me,t1,t2,n1,n2,n3,n4,n12,n14,n24,n23,n13);
+ SPoint2 p1,p2,p3,p4,p13,p24,p23,p14;
reparamMeshEdgeOnFace(n1,n2,gf,p1,p2);
reparamMeshEdgeOnFace(n3,n4,gf,p3,p4);
+ reparamMeshEdgeOnFace(n13,n24,gf,p13,p24);
+ reparamMeshEdgeOnFace(n23,n14,gf,p23,p14);
s_before = surfaceTriangleUV(p1,p2,p4) + surfaceTriangleUV(p1,p2,p3);
s_after = surfaceTriangleUV(p3,p4,p1) + surfaceTriangleUV(p3,p4,p2);
- SPoint2 p34 = (p3+p4)*.5;
- GPoint gp34 = gf->point(p34);
- n34 = new MFaceVertex (gp34.x(),gp34.y(),gp34.z(),gf,p34.x(),p34.y());
+ n34 = optimalLocation (gf,p3,p4);
- const double qold1 = t1->distoShapeMeasure() * t1->gammaShapeMeasure();
- const double qold2 = t2->distoShapeMeasure() * t2->gammaShapeMeasure();
std::vector<MVertex *>vv;
vv.push_back(n13);vv.push_back(n34);vv.push_back(n14);
@@ -314,12 +934,13 @@ struct swap_triangles_p2
vv.push_back(n34);vv.push_back(n23);vv.push_back(n24);
t4 = new MTriangleN (n4,n3,n2,vv,2,t2->getNum(),t2->getPartition());
- const double qnew1 = t3->distoShapeMeasure() * t3->gammaShapeMeasure();
- const double qnew2 = t4->distoShapeMeasure() * t4->gammaShapeMeasure();
+ const double qnew1 = shapeMeasure(t3);
+ const double qnew2 = shapeMeasure(t4);
- quality_old = std::min(qold1,qold2);
- quality_new = std::min(qnew1,qnew2);
+ quality_new = std::min(qnew1,qnew2);
+ quality_old = std::min(qold1,qold2);
+ // printf("old %g linear %g transfinite %g\n",quality_old,quality_new_linear,quality_new_transfinite);
}
bool operator < (const swap_triangles_p2 &other) const
{
@@ -342,6 +963,273 @@ struct swap_triangles_p2
};
+struct swap_triangles_pN
+{
+ MTriangle *t1, *t2;
+ MTriangle *t3, *t4;
+ double quality_old;
+ double quality_new;
+ double s_before,s_after;
+ MVertex *n1, *n2, *n3, *n4;
+ std::vector<MVertex*> n12, n13, n23, n14, n24, n34;
+ std::vector<MVertex*> n123, n124, n134, n234;
+ edgeContainer &edgeVertices;
+ faceContainer &faceVertices;
+ gmshHighOrderSmoother *s;
+
+ swap_triangles_pN(const MEdge &me, MTriangle *_t1, MTriangle *_t2, GFace *gf,
+ edgeContainer &_edgeVertices,
+ faceContainer &_faceVertices,
+ gmshHighOrderSmoother *_s)
+ : t1(_t1), t2(_t2),edgeVertices(_edgeVertices),faceVertices(_faceVertices),s(_s)
+ {
+
+ const double qold1 = shapeMeasure(t1);
+ const double qold2 = shapeMeasure(t2);
+
+ getNodesPN (me,t1,t2,n1,n2,n3,n4,n12,n14,n24,n23,n13);
+ SPoint2 p1,p2,p3,p4,p13,p24,p23,p14;
+ reparamMeshEdgeOnFace(n1,n2,gf,p1,p2);
+ reparamMeshEdgeOnFace(n3,n4,gf,p3,p4);
+
+ s_before = surfaceTriangleUV(p1,p2,p4) + surfaceTriangleUV(p1,p2,p3);
+ s_after = surfaceTriangleUV(p3,p4,p1) + surfaceTriangleUV(p3,p4,p2);
+
+ MTriangle t3lin(n3,n4,n1);
+ MTriangle t4lin(n4,n3,n2);
+
+ t3 = setHighOrder(&t3lin,gf,edgeVertices,faceVertices,false,
+ !t1->getNumFaceVertices(),
+ t1->getPolynomialOrder()-1,s);
+ t4 = setHighOrder(&t4lin,gf,edgeVertices,faceVertices,false,
+ !t1->getNumFaceVertices(),
+ t1->getPolynomialOrder()-1,s);
+
+ optimalLocationPN_ (gf,me, t3, t4,s);
+
+ const double qnew1 = shapeMeasure(t3);
+ const double qnew2 = shapeMeasure(t4);
+
+ quality_new = std::min(qnew1,qnew2);
+ quality_old = std::min(qold1,qold2);
+
+ // if (quality_old < quality_new)
+ printf("QUALITY GOING FROM %12.5E TO %12.5E\n",quality_old,quality_new);
+
+ }
+ bool operator < (const swap_triangles_pN &other) const
+ {
+ return other.quality_new - other.quality_old < quality_new - quality_old;
+ }
+};
+
+
+
+static int optimalLocationP2_ (GFace *gf,
+ const MEdge &me,
+ MTriangle *t1, MTriangle *t2,
+ gmshHighOrderSmoother *s){
+
+ double qini = std::min(shapeMeasure(t1),shapeMeasure(t2));
+
+ if (qini > 0.6) return 0;
+
+ MVertex *n1,*n2,*n3,*n4,*n12,*n14,*n24,*n23,*n13;
+ getNodesP2 (me,t1,t2,n1,n2,n3,n4,n12,n14,n24,n23,n13);
+ SPoint2 p1,p2,p3,p4,p12;
+ reparamMeshVertexOnFace(n12,gf,p12);
+ reparamMeshEdgeOnFace(n3,n4,gf,p3,p4);
+ reparamMeshEdgeOnFace(n3,n4,gf,p3,p4);
+ SPoint2 p34_linear = (p1+p2)*.5;
+ SPoint2 dirt = p4-p3;
+ SPoint2 dirn = (-dirt.y(), dirt.x());
+
+ gmshVector<double> pp(2);
+ pp(0) = p12.x();
+ pp(1) = p12.y();
+ p2data data(gf,t1,t2,n12,s);
+
+ double init = _DeformationEnergy (&data);
+ double opti = minimize_grad_fd (_function_p2tB, pp, &data);
+
+ printf("opti %g init %g\n",opti,init);
+
+ if (init-opti < 1.e-5*(init))return 0;
+
+ double u0 = gf->parBounds(0).low();
+ double u1 = gf->parBounds(0).high();
+ double v0 = gf->parBounds(1).low();
+ double v1 = gf->parBounds(1).high();
+ if (pp(0) < u0 || pp(0) > u1 || pp(1) < v0 || pp(1) > v1)return 0;
+
+ GPoint gp12 = gf->point(SPoint2(pp(0),pp(1)));
+ n12->x() = gp12.x();
+ n12->y() = gp12.y();
+ n12->z() = gp12.z();
+ n12->setParameter(0,pp(0));
+ n12->setParameter(1,pp(1));
+ return 1;
+}
+
+int optimalLocationPN_ (GFace *gf, const MEdge &me, MTriangle *t1, MTriangle *t2,gmshHighOrderSmoother *s){
+ // if quality is sufficient, do nothing (this is an expensive
+ // optimization process)
+ double qopt = std::min(shapeMeasure(t1),shapeMeasure(t2));
+ if (qopt > 0.6) return 0;
+ // get the 2 end vertex of the edge
+ MVertex *n1 = me.getVertex(0);
+ MVertex *n2 = me.getVertex(1);
+ // get all the other nodes that are on the edge
+ int N = t1->getNumVertices();
+ int NF = t1->getNumFaceVertices();
+ int NE = t1->getNumEdgeVertices();
+ std::vector<MVertex*> toOptimize;
+ for (int i=3;i<3+NE;i++){
+ MVertex *v1 = t1->getVertex(i);
+ for (int j=3;j<3+NE;j++){
+ MVertex *v2 = t2->getVertex(j);
+ if (v1 == v2 && v1 != n1 && v1 != n2){
+ toOptimize.push_back(v1);
+ break;
+ }
+ }
+ }
+
+
+ for (int i=3+NE;i<N;i++){
+ toOptimize.push_back(t1->getVertex(i));
+ toOptimize.push_back(t2->getVertex(i));
+ }
+
+ gmshVector<double> pp(2*toOptimize.size());
+ for (int i=0;i<toOptimize.size();i++){
+ SPoint2 pt;
+ reparamMeshVertexOnFace(toOptimize[i],gf,pt);
+ pp(2*i) = pt[0];
+ pp(2*i+1) = pt[1];
+ }
+
+ pNdata data(gf,t1,t2,toOptimize,s);
+ double init = _DeformationEnergy (&data);
+ double opti = minimize_grad_fd (_function_pNtB, pp, &data);
+ ///double opti = minimize_grad_fd (_function_pNt, pp, &data);
+ if (init-opti < 1.e-5*(init))return 0;
+ printf("Optimization has reduced the deformation energy %g -> %g\n",
+ init,opti);
+ for (int i=0;i<toOptimize.size();i++){
+ GPoint gp12 = gf->point(SPoint2(pp(2*i),pp(2*i+1)));
+ toOptimize[i]->x() = gp12.x();
+ toOptimize[i]->y() = gp12.y();
+ toOptimize[i]->z() = gp12.z();
+ toOptimize[i]->setParameter(0,pp(2*i));
+ toOptimize[i]->setParameter(1,pp(2*i+1));
+ }
+ return 1;
+}
+
+
+static int _gmshFindOptimalLocationsP2(GFace *gf, gmshHighOrderSmoother *s)
+{
+ e2t_cont adj;
+ buildEdgeToTriangle(gf->triangles, adj);
+ int N=0;
+ for (e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it){
+ if (it->second.second)
+ N += optimalLocationP2_(gf,it->first, dynamic_cast<MTriangle*>(it->second.first),
+ dynamic_cast<MTriangle*>(it->second.second),s);
+ }
+ return N;
+}
+
+static int _gmshFindOptimalLocationsPN(GFace *gf,gmshHighOrderSmoother *s)
+{
+ printf("coucou1\n");
+ e2t_cont adj;
+ buildEdgeToTriangle(gf->triangles, adj);
+ int N=0;
+ printf("coucou2\n");
+
+ for (e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it){
+ if (it->second.second)
+ N += optimalLocationPN_(gf,it->first, dynamic_cast<MTriangle*>(it->second.first),
+ dynamic_cast<MTriangle*>(it->second.second),s);
+ }
+ printf("coucou3\n");
+ return N;
+}
+
+static int _gmshSwapHighOrderTriangles(GFace *gf,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ gmshHighOrderSmoother *s)
+{
+ e2t_cont adj;
+ buildEdgeToTriangle(gf->triangles, adj);
+
+ std::set<swap_triangles_pN> pairs;
+
+ for (e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it){
+ if (it->second.second){
+ MTriangle *t1 = dynamic_cast<MTriangle*>(it->second.first);
+ MTriangle *t2 = dynamic_cast<MTriangle*>(it->second.second);
+ const double qold1 = shapeMeasure(t1);
+ const double qold2 = shapeMeasure(t2);
+
+ // printf("swap : %g %g\n",qold1,qold2);
+
+ if (qold1 < 0.6 || qold2 < 0.6)
+ pairs.insert(swap_triangles_pN(it->first,t1,t2,gf,
+ edgeVertices,faceVertices,s));
+ }
+ }
+ std::set<swap_triangles_pN>::iterator itp = pairs.begin();
+
+ int nbSwap = 0;
+
+ std::set<MTriangle*> t_removed;
+ std::set<MVertex*> v_removed;
+
+ std::vector<MTriangle*> triangles2;
+ std::vector<MVertex*> mesh_vertices2;
+
+ itp = pairs.begin();
+ while(itp != pairs.end()){
+ double diff = fabs(itp->s_before - itp->s_after);
+ if ( t_removed.find(itp->t1) == t_removed.end() &&
+ t_removed.find(itp->t2) == t_removed.end() &&
+ itp->quality_new > itp->quality_old &&
+ diff < 1.e-9){
+ // itp->print();
+ t_removed.insert(itp->t1);
+ t_removed.insert(itp->t2);
+ triangles2.push_back(itp->t3);
+ triangles2.push_back(itp->t4);
+ // if (itp->n34 != itp->n12){
+ // v_removed.insert(itp->n12);
+ // mesh_vertices2.push_back(itp->n34);
+ // }
+ nbSwap++;
+ }
+ else{
+ delete itp->t3;
+ delete itp->t4;
+ // if (itp->n34 != itp->n12) delete itp->n34;
+ }
+ ++itp;
+ }
+
+ for (unsigned int i = 0; i < gf->triangles.size(); i++){
+ if (t_removed.find(gf->triangles[i]) == t_removed.end()){
+ triangles2.push_back(gf->triangles[i]);
+ }
+ else {
+ delete gf->triangles[i];
+ }
+ }
+ // printf("replacing %d by %d\n",gf->triangles.size(),triangles2.size());
+ gf->triangles = triangles2;
+ return nbSwap;
+}
static int _gmshSwapHighOrderTriangles(GFace *gf)
{
@@ -351,41 +1239,47 @@ static int _gmshSwapHighOrderTriangles(GFace *gf)
std::set<swap_triangles_p2> pairs;
for (e2t_cont::iterator it = adj.begin(); it!= adj.end(); ++it){
- if (it->second.second)
- pairs.insert(swap_triangles_p2(it->first,
- dynamic_cast<MTriangle*>(it->second.first),
- dynamic_cast<MTriangle*>(it->second.second),
- gf));
+ if (it->second.second){
+ MTriangle *t1 = dynamic_cast<MTriangle*>(it->second.first);
+ MTriangle *t2 = dynamic_cast<MTriangle*>(it->second.second);
+ const double qold1 = shapeMeasure(t1);
+ const double qold2 = shapeMeasure(t2);
+ if (qold1 < 0.6 || qold2 < 0.6)
+ pairs.insert(swap_triangles_p2(it->first,t1,t2,gf));
+ }
}
+ std::set<swap_triangles_p2>::iterator itp = pairs.begin();
int nbSwap = 0;
std::set<MTriangle*> t_removed;
std::set<MVertex*> v_removed;
- std::set<swap_triangles_p2>::iterator itp = pairs.begin();
std::vector<MTriangle*> triangles2;
std::vector<MVertex*> mesh_vertices2;
+ itp = pairs.begin();
while(itp != pairs.end()){
double diff = fabs(itp->s_before - itp->s_after);
if ( t_removed.find(itp->t1) == t_removed.end() &&
t_removed.find(itp->t2) == t_removed.end() &&
- itp->quality_new - itp->quality_old > 0.01 &&
+ itp->quality_new > itp->quality_old &&
diff < 1.e-9){
+ // itp->print();
t_removed.insert(itp->t1);
t_removed.insert(itp->t2);
- v_removed.insert(itp->n12);
- itp->print();
triangles2.push_back(itp->t3);
triangles2.push_back(itp->t4);
- mesh_vertices2.push_back(itp->n34);
+ if (itp->n34 != itp->n12){
+ v_removed.insert(itp->n12);
+ mesh_vertices2.push_back(itp->n34);
+ }
nbSwap++;
}
else{
delete itp->t3;
delete itp->t4;
- delete itp->n34;
+ if (itp->n34 != itp->n12) delete itp->n34;
}
++itp;
}
@@ -414,11 +1308,24 @@ static int _gmshSwapHighOrderTriangles(GFace *gf)
return nbSwap;
}
-void gmshSwapHighOrderTriangles(GFace *gf){
- while(_gmshSwapHighOrderTriangles(gf));
+void gmshHighOrderSmoother::swap(GFace *gf,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices){
+ // _gmshSwapHighOrderTriangles(gf);
+ _gmshSwapHighOrderTriangles(gf,edgeVertices,faceVertices,this);
+ //_gmshSwapHighOrderTriangles(gf);
+ //_gmshSwapHighOrderTriangles(gf);
// _gmshSwapHighOrderTriangles(gf);
}
+void gmshHighOrderSmoother::smooth_p2point(GFace *gf){
+ _gmshFindOptimalLocationsP2(gf,this);
+ //_gmshFindOptimalLocationsP2(gf);
+ //_gmshFindOptimalLocationsP2(gf);
+}
+void gmshHighOrderSmoother::smooth_pNpoint(GFace *gf){
+ _gmshFindOptimalLocationsPN(gf,this);
+}
////////////////////////////////////////////////////////////////////////////////
diff --git a/Mesh/gmshSmoothHighOrder.h b/Mesh/gmshSmoothHighOrder.h
index 91858db4b1..9f52604cd3 100644
--- a/Mesh/gmshSmoothHighOrder.h
+++ b/Mesh/gmshSmoothHighOrder.h
@@ -9,27 +9,76 @@
#include <map>
#include <vector>
#include "SVector3.h"
+#include "GmshMatrix.h"
class MVertex;
class MElement;
class GFace;
class GRegion;
-
-void gmshSwapHighOrderTriangles(GFace *gf);
+template<class T> class gmshAssembler;
+class gmshElasticityTerm;
class gmshHighOrderSmoother
{
const int _tag;
- std::map<MVertex*,SVector3> _displ;
+ std::map<MVertex*,SVector3> _straightSidedLocation;
+ std::map<MVertex*,SVector3> _targetLocation;
int _dim;
- void _clean () {_displ.clear();}
+ void _clean () {
+ _straightSidedLocation.clear();
+ _targetLocation.clear();
+ }
+ double _MIDDLE;
+ void moveTo (MVertex *v, const std::map<MVertex*,SVector3> &) const;
public:
gmshHighOrderSmoother (int dim) : _tag(111), _dim(dim) {_clean();}
- void add ( MVertex * v, const SVector3 &d ) {_displ[v] = d;}
+ void add ( MVertex * v, const SVector3 &d ) {
+ _straightSidedLocation[v] = d;
+ _targetLocation[v] = SPoint3(v->x(),v->y(),v->z());
+ }
void smooth ( std::vector<MElement*> & );
- void smooth ( GFace* );
+ double smooth_metric_ ( std::vector<MElement*> &, GFace *gf,
+ gmshAssembler<double> &myAssembler,
+ std::set<MVertex*> &verticesToMove,
+ gmshElasticityTerm &El);
+ void smooth_metric ( std::vector<MElement*> &, GFace *gf );
+ void smooth ( GFace* , bool metric = false);
+ void smooth_p2point ( GFace* );
+ void smooth_pNpoint ( GFace* );
void smooth ( GRegion* );
int getTag() const {return _tag;}
+ void swap(GFace *,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices);
+ void optimize(GFace *,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices);
+ void computeMetricVector (GFace *gf,
+ MElement *e,
+ gmshMatrix<double> &J,
+ gmshMatrix<double> &JT,
+ gmshVector<double> &D);
+ void moveToStraightSidedLocation (MVertex *) const;
+ void moveToTargetLocation (MVertex *) const;
+ void moveToStraightSidedLocation (MElement *) const;
+ void moveToTargetLocation (MElement *) const;
+ void updateTargetLocation (MVertex*, const SPoint3 &, const SPoint2 &) ;
+ inline SVector3 getSSL(MVertex *v) const{
+ std::map<MVertex*,SVector3>::const_iterator it = _straightSidedLocation.find(v);
+ if (it != _straightSidedLocation.end())return it->second;
+ else return SVector3(v->x(),v->y(),v->z());
+ }
+ inline SVector3 getDisplacement(MVertex *v) const{
+ std::map<MVertex*,SVector3>::const_iterator it = _straightSidedLocation.find(v);
+ std::map<MVertex*,SVector3>::const_iterator itt = _targetLocation.find(v);
+ if (it == _straightSidedLocation.end())
+ return SVector3(0.0,0.0,0.0);
+ else{
+ return SVector3(itt->second.x() - it->second.x(),
+ itt->second.y() - it->second.y(),
+ itt->second.z() - it->second.z());
+ }
+ }
};
#endif
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 19cdb668de..55dd9cd302 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -33,30 +33,6 @@ struct xi2lc {
static std::vector<xi2lc> interpLc;
-static void smoothInterpLc(int nbSmooth)
-{
- for(int j = 0; j < nbSmooth; j++){
- for(int i = 0 ; i < (int)interpLc.size(); i++){
- xi2lc &left = (i == 0) ? interpLc[0] : interpLc[i - 1];
- xi2lc &mid = interpLc[i];
- xi2lc &right = (i == (int)interpLc.size() - 1) ?
- interpLc[interpLc.size() - 1] : interpLc[i+1];
- if(1. / mid.lc > 1.1 * 1. / left.lc) mid.lc = left.lc / 1.1;
- if(1. / mid.lc > 1.1 * 1. / right.lc) mid.lc = right.lc / 1.1;
- }
- }
-}
-
-static void printInterpLc(const char *name)
-{
- FILE *f = fopen(name,"w");
- for(unsigned int i = 0; i < interpLc.size(); i++){
- xi2lc &interp = interpLc[i];
- fprintf(f,"%12.5E %12.5E\n", interp.xi, 1 / interp.lc);
- }
- fclose(f);
-}
-
static void buildInterpLc(const std::vector<IntPoint> &lcPoints)
{
IntPoint p;
@@ -74,13 +50,11 @@ static double F_Lc_usingInterpLc(GEdge *ge, double t)
double t1 = it->xi;
double l1 = it->lc;
it++;
- SVector3 der = ge->firstDer(t);
- const double d = norm(der);
- if(it == interpLc.end()) return d * l1;
+ if(it == interpLc.end()) return l1;
double t2 = it->xi;
double l2 = it->lc;
double l = l1 + ((t - t1) / (t2 - t1)) * (l2 - l1);
- return d * l;
+ return l;
}
static double F_Lc_usingInterpLcBis(GEdge *ge, double t)
@@ -92,6 +66,9 @@ static double F_Lc_usingInterpLcBis(GEdge *ge, double t)
double t_begin = bounds.low();
double t_end = bounds.high();
+ SVector3 der = ge->firstDer(t);
+ const double d = norm(der);
+
if(t == t_begin)
lc_here = BGM_MeshSize(ge->getBeginVertex(), t, 0, p.x(), p.y(), p.z());
else if(t == t_end)
@@ -99,7 +76,7 @@ static double F_Lc_usingInterpLcBis(GEdge *ge, double t)
else
lc_here = BGM_MeshSize(ge, t, 0, p.x(), p.y(), p.z());
- return 1 / lc_here;
+ return d / lc_here;
}
static double F_Lc(GEdge *ge, double t)
@@ -305,9 +282,6 @@ void meshGEdge::operator() (GEdge *ge)
Integration(ge, t_begin, t_end, F_Lc_usingInterpLcBis, lcPoints,
CTX::instance()->mesh.lcIntegrationPrecision);
buildInterpLc(lcPoints);
- // printInterpLc("toto1.dat");
- // smoothInterpLc(20);
- // printInterpLc("toto2.dat");
a = Integration(ge, t_begin, t_end, F_Lc_usingInterpLc, Points, 1.e-8);
}
else{
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index de8a3658e6..3c0c4055b6 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -13,6 +13,7 @@
#include "BackgroundMesh.h"
#include "GVertex.h"
#include "GEdge.h"
+#include "GEdgeCompound.h"
#include "GFace.h"
#include "GModel.h"
#include "MVertex.h"
@@ -373,6 +374,26 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER, bool repairSelfInters
v_container all_vertices;
std::map<int, MVertex*>numbered_vertices;
std::list<GEdge*> edges = gf->edges();
+
+ // here, we will replace edges by their compounds
+ if (gf->geomType() == GEntity::CompoundSurface){
+ std::set<GEdge*> mySet;
+ std::list<GEdge*>::iterator it = edges.begin();
+ while(it != edges.end()){
+ if ((*it)->getCompound()){
+ mySet.insert((*it)->getCompound());
+ printf("compound edge %d found in %d\n",(*it)->getCompound()->tag(), (*it)->tag());
+ }
+ else
+ mySet.insert(*it);
+ ++it;
+ }
+ printf("replacing %d edges by %d in the compound %d\n",edges.size(),mySet.size(),gf->tag());
+ edges.clear();
+ edges.insert(edges.begin(), mySet.begin(), mySet.end());
+ }
+
+
std::list<GEdge*> emb_edges = gf->embeddedEdges();
std::list<GVertex*> emb_vertx = gf->embeddedVertices();
std::list<GEdge*>::iterator it = edges.begin();
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index d48963a058..ccac813653 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -32,8 +32,11 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
return l;
}
-inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
- double SCALINGU, double SCALINGV)
+inline double computeEdgeLinearLength(BDS_Point *p1,
+ BDS_Point *p2,
+ GFace *f,
+ double SCALINGU,
+ double SCALINGV)
{
GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
0.5 * (p1->v + p2->v) * SCALINGV));
@@ -49,8 +52,11 @@ inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
return l1 + l2;
}
-inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
- double SCALINGU, double SCALINGV)
+inline double computeEdgeMiddleCoord(BDS_Point *p1,
+ BDS_Point *p2,
+ GFace *f,
+ double SCALINGU,
+ double SCALINGV)
{
if (f->geomType() == GEntity::Plane)
return 0.5;
@@ -73,8 +79,10 @@ inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
return 0.25 * (3 * l2 - l1) / l2;
}
-inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
- double SCALINGU, double SCALINGV)
+inline double computeEdgeLinearLength(BDS_Edge *e,
+ GFace *f,
+ double SCALINGU,
+ double SCALINGV)
{
if (f->geomType() == GEntity::Plane)
return e->length();
@@ -84,7 +92,8 @@ inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
double NewGetLc(BDS_Point *p)
{
- return Extend1dMeshIn2dSurfaces() ? std::min(p->lc(), p->lcBGM()) : p->lcBGM();
+ return Extend1dMeshIn2dSurfaces() ?
+ std::min(p->lc(), p->lcBGM()) : p->lcBGM();
}
double NewGetLc(BDS_Edge *e, GFace *f, double SCALINGU, double SCALINGV)
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 6aebd9d324..9e85a48160 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -14,6 +14,7 @@
#include "meshGFace.h"
#include "GFace.h"
#include "Numeric.h"
+#include "STensor3.h"
const double LIMIT_ = 0.5 * sqrt(2.0);
@@ -521,8 +522,8 @@ static void insertAPoint(GFace *gf,
v->setNum(Us.size());
double lc1 = ((1. - uv[0] - uv[1]) * vSizes[ptin->tri()->getVertex(0)->getNum()] +
- uv[0] * vSizes [ptin->tri()->getVertex(1)->getNum()] +
- uv[1] * vSizes [ptin->tri()->getVertex(2)->getNum()]);
+ uv[0] * vSizes [ptin->tri()->getVertex(1)->getNum()] +
+ uv[1] * vSizes [ptin->tri()->getVertex(2)->getNum()]);
// double eigMetricSurface = gf->getMetricEigenvalue(SPoint2(center[0],center[1]));
double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
vSizesBGM.push_back(lc);
@@ -645,12 +646,29 @@ static double length_metric ( const double p[2], const double q[2], const double
*/
+void testTensor(){
+ SMetric3 t(1.0);
+ t(0,0) = 1;
+ t(1,0) = .2;
+ t(1,1) = 2;
+ t(2,2) = 3;
+ gmshMatrix<double> m(3,3);
+ gmshVector<double> v(3);
+ t.eig(m,v);
+ printf("%12.5E %12.5E %12.5E \n",v(0),v(1),v(2));
+ printf("%12.5E %12.5E %12.5E \n",m(0,0),m(1,0),m(2,0));
+ printf("%12.5E %12.5E %12.5E \n",m(0,1),m(1,1),m(2,1));
+ printf("%12.5E %12.5E %12.5E \n",m(0,2),m(1,2),m(2,2));
+}
+
void gmshBowyerWatsonFrontal(GFace *gf)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
std::set<MTri3*,compareTri3Ptr> ActiveTris;
std::vector<double> vSizes, vSizesBGM, Us, Vs;
+ testTensor();
+
buildMeshGenerationDataStructures(gf, AllTris, vSizes, vSizesBGM, Us, Vs);
// delaunise the initial mesh
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 6a6699c712..92fb0015ba 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -199,20 +199,35 @@ double mesh_functional_distorsion(MTriangle *t, double u, double v)
double normal[3];
prodve(v1b, v2b, normal);
- double sign;
+ double sign = 1.0;
prosca(normal1, normal, &sign);
double det = norm3(normal) * (sign > 0 ? 1. : -1.) / nn;
// compute distorsion
- double dist = std::min(1. / det, det);
- return dist;
+ // double dist = std::min(1. / det, det);
+ return det;
+}
+
+
+double mesh_functional_distorsion_p2(MTriangle *t)
+{
+ double d = mesh_functional_distorsion(t,0.,0.);
+ d = std::min(d,mesh_functional_distorsion(t,1.,0.));
+ d = std::min(d,mesh_functional_distorsion(t,0.,1.));
+ d = std::min(d,mesh_functional_distorsion(t,.5,0.));
+ d = std::min(d,mesh_functional_distorsion(t,0.,0.5));
+ d = std::min(d,mesh_functional_distorsion(t,.5,0.5));
+ return d;
}
+
double qmDistorsionOfMapping (MTriangle *e)
{
// return 1.0;
if (e->getPolynomialOrder() == 1)return 1.0;
+ // if (e->getPolynomialOrder() == 2)return mesh_functional_distorsion_p2(e);
+
IntPt *pts;
int npts;
e->getIntegrationPoints(e->getPolynomialOrder(),&npts, &pts);
@@ -283,5 +298,5 @@ double qmDistorsionOfMapping (MTetrahedron *e)
}
// printf("DMIN = %g\n\n",dmin);
- return dmin< 0 ? 0 :dmin;
+ return dmin;
}
diff --git a/Numeric/GmshMatrix.cpp b/Numeric/GmshMatrix.cpp
index b43788bf68..96281ed051 100644
--- a/Numeric/GmshMatrix.cpp
+++ b/Numeric/GmshMatrix.cpp
@@ -104,8 +104,68 @@ extern "C" {
void dgesvd_(const char* jobu, const char *jobvt, int *M, int *N,
double *A, int *lda, double *S, double* U, int *ldu,
double *VT, int *ldvt, double *work, int *lwork, int *info);
+ void dgeev_(const char *jobvl, const char *jobvr,
+ int *n, double *a, int *lda,
+ double *wr, double *wi,
+ double *vl, int *ldvl,
+ double *vr, int *ldvr,
+ double *work, int *lwork,
+ int *info);
}
+template<>
+bool gmshMatrix<double>::invertInPlace()
+{
+ int N = size1(), nrhs = N, lda = N, ldb = N, info;
+ int *ipiv = new int[N];
+ double * invA = new double[N*N];
+
+ for (size_t i=0;i<N*N;i++) invA[i ] = 0.;
+ for (size_t i=0;i<N;i++) invA[i*N+i] = 1.;
+
+ dgesv_(&N, &nrhs, _data, &lda, ipiv, invA, &ldb, &info);
+ std::memcpy(_data,invA,N*N*sizeof(double));
+
+ delete [] invA;
+ delete [] ipiv;
+
+ if(info == 0) return true;
+ if(info > 0)
+ Msg::Error("U(%d,%d)=0 in matric inversion", info, info);
+ else
+ Msg::Error("Wrong %d-th argument in matrix inversion", -info);
+ return false;
+}
+
+
+template<>
+bool gmshMatrix<double>::eig(gmshMatrix<double> &VL, // left eigenvectors
+ gmshVector<double> &DR, // Real part of eigenvalues
+ gmshVector<double> &DI, // Im part of eigenvalues
+ gmshMatrix<double> &VR )
+{
+ int N = size1(), info;
+ int LWORK = 10*N;
+ double * work = new double[LWORK];
+
+ dgeev_("V","V",
+ &N,_data,
+ &N,DR._data,DI._data,
+ VL._data,&N,
+ VR._data,&N,
+ work,&LWORK,&info);
+
+ delete [] work;
+
+ if(info == 0) return true;
+ if(info > 0)
+ Msg::Error("QR Algorithm failed to compute all the eigenvalues", info, info);
+ else
+ Msg::Error("Wrong %d-th argument in eig", -info);
+ return false;
+}
+
+
template<>
bool gmshMatrix<double>::lu_solve(const gmshVector<double> &rhs, gmshVector<double> &result)
{
diff --git a/Numeric/GmshMatrix.h b/Numeric/GmshMatrix.h
index 40df2676e8..771f963b79 100644
--- a/Numeric/GmshMatrix.h
+++ b/Numeric/GmshMatrix.h
@@ -53,6 +53,14 @@ class gmshVector
else
for(int i = 0; i < _r; ++i) _data[i] *= s;
}
+
+ inline scalar operator *(const gmshVector<scalar> & other)
+ {
+ scalar s = 0.0;
+ for(int i = 0; i < _r; ++i) s += _data[i]*other._data[i];
+ return s;
+ }
+
};
template <class scalar>
@@ -173,6 +181,25 @@ class gmshMatrix
Msg::Error("LU factorization requires LAPACK");
return false;
}
+#endif
+ ;
+ bool invertInPlace()
+#if !defined(HAVE_LAPACK)
+ {
+ Msg::Error("Matrix inversion requires LAPACK");
+ return false;
+ }
+#endif
+ ;
+ bool eig(gmshMatrix<scalar> &VL, // left eigenvectors
+ gmshVector<double> &DR, // Real part of eigenvalues
+ gmshVector<double> &DI, // Im part of eigen
+ gmshMatrix<scalar> &VR)
+#if !defined(HAVE_LAPACK)
+ {
+ Msg::Error("Eigenvalue computations requires LAPACK");
+ return false;
+ }
#endif
;
bool invert(gmshMatrix<scalar> &result)
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 6d6804c3e0..8519c76896 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -571,3 +571,125 @@ bool newton_fd(void (*func)(gmshVector<double> &, gmshVector<double> &, void *),
}
return false;
}
+
+/*
+min_a f(x+a*d);
+
+f(x+a*d) = f(x) + f'(x) (
+
+*/
+
+void gmshLineSearch (double (*func)(gmshVector<double> &, void *), void* data,
+ gmshVector<double> &x,
+ gmshVector<double> &p,
+ gmshVector<double> &g,
+ double &f,
+ double stpmax, int &check)
+{
+ int i;
+ double alam,alam2,alamin,f2,fold2,rhs1,rhs2,sum,temp,
+ tmplam;
+
+ const double ALF = 1.0e-4;
+ const double TOLX = 1.0e-9;
+
+ gmshVector<double> xold(x);
+ const double fold = (*func)(xold,data);
+
+ check=0;
+ int n = x.size();
+ double norm = p.norm();
+ if (norm > stpmax)p.scale(stpmax/norm);
+ double slope=0.0;
+ for (i=0;i<n;i++)slope += g(i)*p(i);
+ double test=0.0;
+ for (i=0;i<n;i++) {
+ temp=fabs(p(i))/std::max(fabs(xold(i)),1.0);
+ if (temp > test) test=temp;
+ }
+ /*
+ for (int j=0;j<100;j++){
+ double sx = (double)j/99;
+ for (i=0;i<n;i++) x(i)=xold(i)+10*sx*p(i);
+ double jzede = (*func)(x,data);
+ }
+ */
+
+ alamin=TOLX/test;
+ alam=1.0;
+ while(1) {
+ for (i=0;i<n;i++) x(i)=xold(i)+alam*p(i);
+ f=(*func)(x,data);
+ // printf("f = %g x = %g %g alam = %g p = %g %g\n",f,x(0),x(1),alam,p(0),p(1));
+ if (alam < alamin) {
+ for (i=0;i<n;i++) x(i)=xold(i);
+ // printf("ALERT : alam %g alamin %g\n",alam,alamin);
+ check=1;
+ return;
+ }
+ else if (f <= fold + ALF*alam*slope) return;
+ else {
+ if (alam == 1.0)
+ tmplam = -slope/(2.0*(f-fold-slope));
+ else {
+ rhs1 = f-fold-alam*slope;
+ rhs2=f2-fold2-alam2*slope;
+ const double a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
+ const double b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
+ if (a == 0.0) tmplam = -slope/(2.0*b);
+ else {
+ const double disc=b*b-3.0*a*slope;
+ if (disc<0.0) Msg::Error("Roundoff problem in gmshLineSearch.");
+ else tmplam=(-b+sqrt(disc))/(3.0*a);
+ }
+ if (tmplam>0.5*alam)
+ tmplam=0.5*alam;
+ }
+ }
+ alam2=alam;
+ f2 = f;
+ fold2=fold;
+ alam=std::max(tmplam,0.1*alam);
+ }
+}
+
+double minimize_grad_fd (double (*func)(gmshVector<double> &, void *),
+ gmshVector<double> &x, void *data)
+{
+ const int MAXIT = 3;
+ const double EPS = 1.e-4;
+ const int N = x.size();
+
+ gmshVector<double> grad(N);
+ gmshVector<double> dir(N);
+ double f,feps,finit;
+
+ for (int iter = 0; iter < MAXIT; iter++){
+ // compute gradient of func
+ f = func(x,data);
+ if (iter == 0)finit = f;
+ // printf("Opti iter %d x = (%g %g) f = %g\n",iter,x(0),x(1),f);
+ // printf("grad = (");
+ for (int j = 0; j < N; j++){
+ double h = EPS * fabs(x(j));
+ if(h == 0.) h = EPS;
+ x(j) += h;
+ feps = func(x, data);
+ grad(j) = (feps - f) / h;
+ // printf("%g ",grad(j));
+ dir(j) = -grad(j);
+ x(j) -= h;
+ }
+ // printf(")\n ");
+ // do a 1D line search to fine the minimum
+ // of f(x - \alpha \nabla f)
+ double f, stpmax=100000;
+ int check;
+ gmshLineSearch (func, data, x,dir,grad,f,stpmax,check);
+ // printf("Line search done x = (%g %g) f = %g\n",x(0),x(1),f);
+ if (check == 1)break;
+ }
+
+ return f;
+}
+
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 578a66f736..ba4e0deb72 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -74,5 +74,7 @@ double ComputeScalarRep(int numComp, double *val);
void invert_singular_matrix3x3(double MM[3][3], double II[3][3]);
bool newton_fd(void (*func)(gmshVector<double> &, gmshVector<double> &, void *),
gmshVector<double> &x, void *data, double relax=1., double tolx=1.e-6);
+double minimize_grad_fd (double (*func)(gmshVector<double> &, void *),
+ gmshVector<double> &x, void *data);
#endif
diff --git a/Numeric/gmshAssembler.h b/Numeric/gmshAssembler.h
index bad4d5b052..56d56d9be7 100644
--- a/Numeric/gmshAssembler.h
+++ b/Numeric/gmshAssembler.h
@@ -166,6 +166,13 @@ class gmshAssembler {
}
int sizeOfR() const { return numbering.size(); }
int sizeOfF() const { return fixed.size(); }
+ void systemSolve(){
+ lsys->systemSolve();
+ }
+ void systemClear(){
+ lsys->zeroMatrix();
+ lsys->zeroRightHandSide();
+ }
};
#endif
diff --git a/Numeric/gmshElasticity.cpp b/Numeric/gmshElasticity.cpp
index 2750a47ff1..3d323ca272 100644
--- a/Numeric/gmshElasticity.cpp
+++ b/Numeric/gmshElasticity.cpp
@@ -9,7 +9,7 @@
void gmshElasticityTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
{
int nbNodes = e->getNumVertices();
- int integrationOrder = 3 * (e->getPolynomialOrder() - 1);
+ int integrationOrder = 3 * (e->getPolynomialOrder() - 1) ;
int npts;
IntPt *GP;
double jac[3][3];
diff --git a/Numeric/gmshElasticity.h b/Numeric/gmshElasticity.h
index 037aeedbd2..ac4c81268c 100644
--- a/Numeric/gmshElasticity.h
+++ b/Numeric/gmshElasticity.h
@@ -13,10 +13,11 @@
#include "GmshMatrix.h"
class gmshElasticityTerm : public gmshNodalFemTerm<double> {
+ protected:
double _E, _nu;
int _iField;
SVector3 _f;
- protected:
+ public:
virtual int sizeOfR(MElement *e) const { return 3 * e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return 3 * e->getNumVertices(); }
void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
diff --git a/Numeric/gmshLaplace.cpp b/Numeric/gmshLaplace.cpp
index 0c9f54af54..d368962cd5 100644
--- a/Numeric/gmshLaplace.cpp
+++ b/Numeric/gmshLaplace.cpp
@@ -50,3 +50,7 @@ void gmshLaplaceTerm::elementMatrix(MElement *e, gmshMatrix<double> &m) const
for (int k = 0; k < j; k++)
m(k, j) = m(j, k);
}
+
+void gmshLaplaceTerm2DParametric::elementMatrix(MElement *e, gmshMatrix<double> &m) const
+{
+}
diff --git a/Numeric/gmshLaplace.h b/Numeric/gmshLaplace.h
index 3a462294f9..223b5696da 100644
--- a/Numeric/gmshLaplace.h
+++ b/Numeric/gmshLaplace.h
@@ -14,10 +14,10 @@
#include "GmshMatrix.h"
class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
- private:
+ protected:
const gmshFunction<double> *_diffusivity;
const int _iField ;
- protected:
+ public:
virtual int sizeOfR(MElement *e) const { return e->getNumVertices(); }
virtual int sizeOfC(MElement *e) const { return e->getNumVertices(); }
void getLocalDofR(MElement *e, int iRow, MVertex **vR, int *iCompR, int *iFieldR) const
@@ -32,4 +32,13 @@ class gmshLaplaceTerm : public gmshNodalFemTerm<double> {
virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
};
+class gmshLaplaceTerm2DParametric : gmshLaplaceTerm {
+ GFace *_gf;
+ public:
+ gmshLaplaceTerm2DParametric(GFace *gf, gmshFunction<double> *diffusivity, int iField = 0) :
+ _gf(gf),gmshLaplaceTerm(gf->model(),diffusivity,iField){}
+ virtual void elementMatrix(MElement *e, gmshMatrix<double> &m) const;
+};
+
+
#endif
diff --git a/Numeric/gmshLinearSystemGmm.h b/Numeric/gmshLinearSystemGmm.h
index a1caffafd0..1ae20d1e84 100644
--- a/Numeric/gmshLinearSystemGmm.h
+++ b/Numeric/gmshLinearSystemGmm.h
@@ -75,8 +75,8 @@ class gmshLinearSystemGmm : public gmshLinearSystem<scalar> {
void setGmres(int n){ _gmres = n; }
virtual int systemSolve()
{
- //gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 10, 0.);
- gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 2, 1.e-10);
+ //gmm::ilutp_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 15, 0.);
+ gmm::ildltt_precond<gmm::row_matrix<gmm::wsvector<scalar> > > P(*_a, 10, 1.e-10);
gmm::iteration iter(_prec);
iter.set_noisy(_noisy);
if(_gmres) gmm::gmres(*_a, *_x, *_b, P, 100, iter);
diff --git a/Numeric/gmshTermOfFormulation.h b/Numeric/gmshTermOfFormulation.h
index 2ee14a8d3c..d42d47817b 100644
--- a/Numeric/gmshTermOfFormulation.h
+++ b/Numeric/gmshTermOfFormulation.h
@@ -29,7 +29,7 @@ class gmshTermOfFormulation {
// of the mesh
template<class scalar>
class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
- protected:
+ public:
// return the number of columns of the element matrix
virtual int sizeOfC(MElement*) const = 0;
// return the number of rows of the element matrix
@@ -181,7 +181,6 @@ class gmshNodalFemTerm : public gmshTermOfFormulation<scalar> {
}
}
}
-
};
#endif
diff --git a/Parser/Makefile b/Parser/Makefile
index 446dfe2d54..d29a736f87 100644
--- a/Parser/Makefile
+++ b/Parser/Makefile
@@ -73,9 +73,9 @@ Gmsh.tab${OBJEXT}: Gmsh.tab.cpp ../Common/GmshConfig.h ../Common/GmshMessage.h \
../Fltk/Draw.h ../Common/Options.h ../Post/ColorTable.h \
../Common/Colors.h ../Common/Options.h Parser.h ../Common/OpenFile.h \
../Common/CommandLine.h FunctionManager.h ../Common/OS.h \
- ../Common/CreateFile.h ../Mesh/Field.h ../Post/PView.h \
- ../Mesh/BackgroundMesh.h ../Post/PViewDataList.h ../Post/PViewData.h \
- ../Plugin/PluginManager.h
+ ../Common/CreateFile.h ../Mesh/Field.h ../Geo/STensor3.h \
+ ../Geo/SVector3.h ../Post/PView.h ../Mesh/BackgroundMesh.h \
+ ../Post/PViewDataList.h ../Post/PViewData.h ../Plugin/PluginManager.h
Gmsh.yy${OBJEXT}: Gmsh.yy.cpp ../Common/GmshMessage.h ../Geo/Geo.h \
../Common/GmshDefines.h ../Geo/gmshSurface.h ../Geo/Pair.h \
../Geo/Range.h ../Geo/SPoint2.h ../Geo/SPoint3.h ../Geo/SVector3.h \
diff --git a/Plugin/Makefile b/Plugin/Makefile
index 332c8d5a6a..5e6110d294 100644
--- a/Plugin/Makefile
+++ b/Plugin/Makefile
@@ -199,18 +199,20 @@ Skin${OBJEXT}: Skin.cpp Skin.h Plugin.h ../Common/Options.h ../Post/ColorTable.h
../Common/ListUtils.h ../Common/TreeUtils.h ../Common/avl.h \
../Common/ListUtils.h ../Common/MallocUtils.h ../Common/Context.h \
../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
-GSHHS${OBJEXT}: GSHHS.cpp ../Mesh/Field.h ../Common/GmshConfig.h ../Post/PView.h \
- ../Geo/SPoint3.h ../Geo/SPoint2.h GSHHS.h Plugin.h ../Common/Options.h \
- ../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PViewDataList.h \
- ../Post/PViewData.h ../Geo/SBoundingBox3d.h ../Geo/SPoint3.h \
- ../Numeric/GmshMatrix.h ../Common/ListUtils.h ../Geo/GModel.h \
- ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+GSHHS${OBJEXT}: GSHHS.cpp ../Mesh/Field.h ../Common/GmshConfig.h \
+ ../Geo/STensor3.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+ ../Numeric/GmshMatrix.h ../Common/GmshMessage.h ../Numeric/Numeric.h \
+ ../Numeric/GmshMatrix.h ../Post/PView.h ../Geo/SPoint2.h GSHHS.h \
+ Plugin.h ../Common/Options.h ../Post/ColorTable.h \
+ ../Post/PViewDataList.h ../Post/PViewData.h ../Geo/SBoundingBox3d.h \
+ ../Geo/SPoint3.h ../Common/ListUtils.h ../Geo/GModel.h ../Geo/GVertex.h \
+ ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
../Geo/SBoundingBox3d.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/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
Extract${OBJEXT}: Extract.cpp ../Common/GmshConfig.h Extract.h Plugin.h \
../Common/Options.h ../Post/ColorTable.h ../Common/GmshMessage.h \
../Post/PView.h ../Geo/SPoint3.h ../Post/PViewDataList.h \
@@ -244,13 +246,14 @@ FieldView${OBJEXT}: FieldView.cpp FieldView.h Plugin.h ../Common/Options.h \
../Geo/SPoint3.h ../Post/PViewDataList.h ../Post/PViewData.h \
../Geo/SBoundingBox3d.h ../Geo/SPoint3.h ../Numeric/GmshMatrix.h \
../Common/GmshConfig.h ../Common/ListUtils.h ../Mesh/Field.h \
- ../Geo/GModel.h ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h \
- ../Geo/SPoint3.h ../Geo/SBoundingBox3d.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/STensor3.h ../Geo/SVector3.h ../Geo/SPoint3.h \
+ ../Numeric/Numeric.h ../Numeric/GmshMatrix.h ../Geo/GModel.h \
+ ../Geo/GVertex.h ../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
+ ../Geo/SBoundingBox3d.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
Integrate${OBJEXT}: Integrate.cpp Integrate.h Plugin.h ../Common/Options.h \
../Post/ColorTable.h ../Common/GmshMessage.h ../Post/PView.h \
diff --git a/Post/shapeFunctions.h b/Post/shapeFunctions.h
index 965347251d..1871bfdc16 100644
--- a/Post/shapeFunctions.h
+++ b/Post/shapeFunctions.h
@@ -84,6 +84,7 @@ public:
}
return sqrt(SQU(jac[0][0])+SQU(jac[0][1])+SQU(jac[0][2]));
default:
+ jac[0][0] = jac[1][1] = jac[2][2] = 1.;
return 1.;
}
}
diff --git a/benchmarks/2d/naca12_2d.geo b/benchmarks/2d/naca12_2d.geo
index 3b2bd214b3..6baf8b1e2f 100644
--- a/benchmarks/2d/naca12_2d.geo
+++ b/benchmarks/2d/naca12_2d.geo
@@ -1,6 +1,6 @@
-lc = 0.2;
-lc2 = 1;
-lc3 = 0.1;
+lc = 0.2 ;
+lc2 = 1 ;
+lc3 = 0.1 ;
Point(1) = {1.000000e+00,0.000000e+00,0.000000e+00,lc3};
Point(2) = {9.997533e-01,0.000000e+00,-3.498543e-05,lc};
Point(3) = {9.990134e-01,0.000000e+00,-1.398841e-04,lc};
@@ -224,7 +224,17 @@ Line Loop(9) = {6,7,8,5};
Line Loop(10) = {2,3,4,1};
Plane Surface(11) = {9,10};
+Physical Surface(11)={11};
//Point(9999) = {0.6,0,0,1};
-//Recombine Surface {11};
-
+Field[1] = Attractor;
+Field[1].EdgesList = {1, 2, 3, 4};
+Field[1].NNodesByEdge = 2000;
+Field[2] = Laplacian;
+Field[2].Delta = 0.001;
+Field[2].IField = 1;
+Field[3] = MathEval;
+Field[3].F = "20/(F2*3.14)";
+//Background Field = 3;
+Field[4] = Gradient;
+Field[4].Delta = 0.001;
--
GitLab