diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index cc556d3d120dfc465907d84f4619b5b994405c8a..c3348d7ec584abf14ec9c7cced1c9820ff275eee 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -27,6 +27,7 @@
#if defined(HAVE_OPTHOM)
#include "OptHomRun.h"
+#include "OptHomElastic.h"
#endif
static void change_completeness_cb(Fl_Widget *w, void *data)
@@ -115,11 +116,11 @@ static void highordertools_runopti_cb(Fl_Widget *w, void *data)
int nbLayers = (int) o->value[2]->value();
double threshold_max = o->value[8]->value();
+#if defined(HAVE_OPTHOM)
if(elastic){
ElasticAnalogy(GModel::current(), threshold_min, onlyVisible);
}
- else {
-#if defined(HAVE_OPTHOM)
+ else{
OptHomParameters p;
p.nbLayers = nbLayers;
p.BARRIER_MIN = threshold_min;
@@ -137,8 +138,10 @@ static void highordertools_runopti_cb(Fl_Widget *w, void *data)
p.adaptBlobLayerFact = (int) o->value[10]->value();
p.adaptBlobDistFact = o->value[11]->value();
HighOrderMeshOptimizer (GModel::current(),p);
-#endif
}
+#else
+ Msg::Error("High-order mesh optimization requires the OPTHOM module");
+#endif
CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE | ENT_VOLUME);
drawContext::global()->draw();
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 5d365bec22f221d128ec4abf5913d79bcc87072f..fd4bc45550c8d07bec4fc5169ba8aa6f71800d7d 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -23,7 +23,6 @@ set(SRC
BoundaryLayers.cpp
BDS.cpp
HighOrder.cpp
- highOrderTools.cpp
meshPartition.cpp
meshRefine.cpp
multiscalePartition.cpp
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index ffcf0b3b23cbdfbeb55c22e93cef910fdb0d8854..f75ee9cb8da46a38a785b7f99faa4e6eb2f89f1a 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -8,7 +8,6 @@
//
#include "HighOrder.h"
-#include "highOrderTools.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
@@ -100,7 +99,6 @@ static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
return false;
}
-
static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
double v0, double vN, int N,
double *u, double *v)
@@ -121,7 +119,6 @@ static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
t[N] = 1.0;
return true;
-
}
static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
@@ -685,10 +682,10 @@ 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)
+static MTriangle *setHighOrder(MTriangle *t, GFace *gf,
+ edgeContainer &edgeVertices,
+ faceContainer &faceVertices,
+ bool linear, bool incomplete, int nPts)
{
std::vector<MVertex*> ve, vf;
getEdgeVertices(gf, t, ve, edgeVertices, linear, nPts);
@@ -898,36 +895,28 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
}
gr->prisms = prisms2;
-/* * * * * * * * * * * * * * * * * PYRAMIDS * * * * * * * * * * * * * * * * * */
-
std::vector<MPyramid*> pyramids2;
-
for(unsigned int i = 0; i < gr->pyramids.size(); i++) {
MPyramid *p = gr->pyramids[i];
std::vector<MVertex*> ve, vf, vr;
getEdgeVertices(gr, p, ve, edgeVertices, linear, nPts);
-
getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
-
ve.insert(ve.end(), vf.begin(), vf.end());
-
vr.reserve((nPts-1)*(nPts)*(2*(nPts-1)+1)/6);
-
int verts_lvl3[12] = {37,40,38,43,46,44,49,52,50,55,58,56};
-
int verts_lvl2[8];
if (nPts == 4) {
verts_lvl2[0] = 42; verts_lvl2[1] = 41;
verts_lvl2[2] = 48; verts_lvl2[3] = 47;
verts_lvl2[4] = 54; verts_lvl2[5] = 53;
verts_lvl2[6] = 60; verts_lvl2[7] = 59;
- } else {
+ }
+ else {
verts_lvl2[0] = 29; verts_lvl2[1] = 30;
verts_lvl2[2] = 35; verts_lvl2[3] = 36;
verts_lvl2[4] = 38; verts_lvl2[5] = 39;
verts_lvl2[6] = 32; verts_lvl2[7] = 33;
}
-
int verts_lvl1[4];
switch(nPts) {
case(4):
@@ -949,7 +938,6 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
verts_lvl1[3] = 22;
break;
}
-
for (int q = 0; q < nPts - 1; q++) {
std::vector<MVertex*> vq, veq;
vq.push_back(ve[2*nPts + q]);
@@ -980,16 +968,15 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
vr.insert(cursor, v);
}
else if (nPts-q == 3) {
-
MQuadrangleN incpl2(vq[0], vq[1], vq[2], vq[3], veq, 3);
int offsets[4] = {nPts == 4 ? 7 : 0,
nPts == 4 ? 9 : 1,
nPts == 4 ? 11 : 2,
nPts == 4 ? 12 : 3};
double quad_v [4][2] = {{-1.0/3.0, -1.0/3.0},
- { 1.0/3.0, -1.0/3.0},
- { 1.0/3.0, 1.0/3.0},
- {-1.0/3.0, 1.0/3.0}};
+ { 1.0/3.0, -1.0/3.0},
+ { 1.0/3.0, 1.0/3.0},
+ {-1.0/3.0, 1.0/3.0}};
SPoint3 pointz;
for (int k = 0; k<4; k++) {
incpl2.pnt(quad_v[k][0], quad_v[k][1], 0, pointz);
@@ -1004,16 +991,16 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
MQuadrangleN incpl2(vq[0], vq[1], vq[2], vq[3], veq, 4);
int offsets[9] = {0, 1, 2, 3, 5, 8, 10, 6, 13};
double quad_v [9][2] = {
- { -0.5, -0.5},
- { 0.5, -0.5},
- { 0.5, 0.5},
- { -0.5, 0.5},
- { 0.0, -0.5},
- { 0.5, 0.0},
- { 0.0, 0.5},
- { -0.5, 0.0},
- { 0.0, 0.0}
- };
+ { -0.5, -0.5},
+ { 0.5, -0.5},
+ { 0.5, 0.5},
+ { -0.5, 0.5},
+ { 0.0, -0.5},
+ { 0.5, 0.0},
+ { 0.0, 0.5},
+ { -0.5, 0.0},
+ { 0.0, 0.0}
+ };
SPoint3 pointz;
for (int k = 0; k<9; k++) {
incpl2.pnt(quad_v[k][0], quad_v[k][1], 0, pointz);
@@ -1024,25 +1011,21 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
vr.insert(cursor, v);
}
}
-
}
ve.insert(ve.end(), vr.begin(), vr.end());
MPyramid *n = new MPyramidN(p->getVertex(0), p->getVertex(1),
p->getVertex(2), p->getVertex(3),
p->getVertex(4), ve, nPts + 1,
- 0, p->getPartition());
+ 0, p->getPartition());
pyramids2.push_back(n);
SPoint3 test_pnt;
n->pnt(-1,-1,0, test_pnt);
-
delete p;
}
gr->pyramids = pyramids2;
gr->deleteVertexArrays();
}
-/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
-
template<class T>
static void setFirstOrder(GEntity *e, std::vector<T*> &elements, bool onlyVisible)
{
@@ -1148,8 +1131,8 @@ void checkHighOrderTriangles(const char* cc, GModel *m,
minGGlob, avg / (count ? count : 1));
}
-static void checkHighOrderTetrahedron(const char* cc, GModel *m,
- std::vector<MElement*> &bad, double &minJGlob)
+void checkHighOrderTetrahedron(const char* cc, GModel *m,
+ std::vector<MElement*> &bad, double &minJGlob)
{
bad.clear();
minJGlob = 1.0;
@@ -1176,58 +1159,6 @@ static void checkHighOrderTetrahedron(const char* cc, GModel *m,
avg / (count ? count : 1));
}
-extern double mesh_functional_distorsion_2D(MElement *t, double u, double v);
-
-void printJacobians(GModel *m, const char *nm)
-{
- const int n = 100;
- double D[n][n], X[n][n], Y[n][n], Z[n][n];
-
- FILE *f = Fopen(nm,"w");
- fprintf(f,"View \"\"{\n");
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
- for(unsigned int j = 0; j < (*it)->triangles.size(); j++){
- MTriangle *t = (*it)->triangles[j];
- for(int i = 0; i < n; i++){
- for(int k = 0; k < n - i; k++){
- SPoint3 pt;
- double u = (double)i / (n - 1);
- double v = (double)k / (n - 1);
- t->pnt(u, v, 0, pt);
- D[i][k] = 0.; //mesh_functional_distorsion_2D(t, u, v);
- //X[i][k] = u;
- //Y[i][k] = v;
- //Z[i][k] = 0.0;
- X[i][k] = pt.x();
- Y[i][k] = pt.y();
- Z[i][k] = pt.z();
- }
- }
- for(int i= 0; i < n -1; i++){
- for(int k = 0; k < n - i -1; k++){
- fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
- X[i][k],Y[i][k],Z[i][k],
- X[i+1][k],Y[i+1][k],Z[i+1][k],
- X[i][k+1],Y[i][k+1],Z[i][k+1],
- D[i][k],
- D[i+1][k],
- D[i][k+1]);
- if (i != n-2 && k != n - i -2)
- fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
- X[i+1][k],Y[i+1][k],Z[i+1][k],
- X[i+1][k+1],Y[i+1][k+1],Z[i+1][k+1],
- X[i][k+1],Y[i][k+1],Z[i][k+1],
- D[i+1][k],
- D[i+1][k+1],
- D[i][k+1]);
- }
- }
- }
- }
- fprintf(f,"};\n");
- fclose(f);
-}
-
void getMeshInfoForHighOrder(GModel *gm, int &meshOrder, bool &complete,
bool &CAD)
{
@@ -1252,44 +1183,6 @@ void getMeshInfoForHighOrder(GModel *gm, int &meshOrder, bool &complete,
}
}
-void ElasticAnalogy(GModel *m, double threshold, bool onlyVisible)
-{
- bool CAD, complete;
- int meshOrder;
-
- getMeshInfoForHighOrder (m,meshOrder, complete, CAD);
- highOrderTools hot(m);
- // now we smooth mesh the internal vertices of the faces
- // we do that model face by model face
- std::vector<MElement*> bad;
- double worst;
- checkHighOrderTriangles("Surface mesh", m, bad, worst);
- {
- for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
- if (onlyVisible && !(*it)->getVisibility()) continue;
- std::vector<MElement*> v;
- v.insert(v.begin(), (*it)->triangles.begin(), (*it)->triangles.end());
- v.insert(v.end(), (*it)->quadrangles.begin(), (*it)->quadrangles.end());
- if (CAD) hot.applySmoothingTo(v, (*it));
- else hot.applySmoothingTo(v, 1.e32, false);
- }
- }
- checkHighOrderTriangles("Final surface mesh", m, bad, worst);
-
- checkHighOrderTetrahedron("Volume Mesh", m, bad, worst);
- {
- for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
- if (onlyVisible && !(*it)->getVisibility())continue;
- std::vector<MElement*> v;
- v.insert(v.begin(), (*it)->tetrahedra.begin(), (*it)->tetrahedra.end());
- v.insert(v.end(), (*it)->hexahedra.begin(), (*it)->hexahedra.end());
- v.insert(v.end(), (*it)->prisms.begin(), (*it)->prisms.end());
- hot.applySmoothingTo(v,1.e32,false);
- }
- }
- checkHighOrderTetrahedron("File volume Mesh", m, bad, worst);
-}
-
void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisible)
{
// replace all the elements in the mesh with second order elements
@@ -1364,7 +1257,7 @@ void computeDistanceFromMeshToGeometry (GModel *m, distanceFromMeshToGeometry_t
{
for (GModel::eiter itEdge = m->firstEdge(); itEdge != m->lastEdge(); ++itEdge) {
double d2,dmax;
- (*itEdge)->computeDistanceFromMeshToGeometry (d2,dmax);
+ (*itEdge)->computeDistanceFromMeshToGeometry(d2, dmax);
dist.d2[*itEdge] = d2;
dist.d_max[*itEdge] = dmax;
}
diff --git a/Mesh/HighOrder.h b/Mesh/HighOrder.h
index f7281a54610c861117ce5cf376b63ed8c8891540..b4641e81f1611b7cd6461ca952a68e009f5f2e2c 100644
--- a/Mesh/HighOrder.h
+++ b/Mesh/HighOrder.h
@@ -22,15 +22,14 @@ typedef std::map<MFace, std::vector<MVertex*>, Less_Face> faceContainer;
void SetOrder1(GModel *m, bool onlyVisible = false);
void SetOrderN(GModel *m, int order, bool linear=true, bool incomplete=false,
bool onlyVisible=false);
-void ElasticAnalogy(GModel *m, double threshold, bool onlyVisible);
-void SetHighOrderComplete (GModel *m, bool onlyVisible);
-void SetHighOrderInComplete (GModel *m, bool onlyVisible);
-MTriangle* setHighOrder(MTriangle *t, GFace *gf,
- edgeContainer &edgeVertices,
- faceContainer &faceVertices,
- bool linear, bool incomplete, int nPts = 1);
+
+void SetHighOrderComplete(GModel *m, bool onlyVisible);
+void SetHighOrderInComplete(GModel *m, bool onlyVisible);
+
void checkHighOrderTriangles(const char* cc, GModel *m,
std::vector<MElement*> &bad, double &minJGlob);
+void checkHighOrderTetrahedron(const char* cc, GModel *m,
+ std::vector<MElement*> &bad, double &minJGlob);
struct distanceFromMeshToGeometry_t {
std::map<GEntity*, double> d_max, d2;
diff --git a/contrib/HighOrderMeshOptimizer/CMakeLists.txt b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
index 6494efaef68843971eae56fc9bf19c1fd2f7037d..ffe20b2110c3cd6700608355b48ceeced33d071a 100644
--- a/contrib/HighOrderMeshOptimizer/CMakeLists.txt
+++ b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
@@ -8,6 +8,7 @@ set(SRC
OptHOM.cpp
OptHomRun.cpp
ParamCoord.cpp
+ OptHomElastic.cpp
)
file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.hpp)
diff --git a/Mesh/highOrderTools.cpp b/contrib/HighOrderMeshOptimizer/OptHomElastic.cpp
similarity index 87%
rename from Mesh/highOrderTools.cpp
rename to contrib/HighOrderMeshOptimizer/OptHomElastic.cpp
index 0ff338032e0705ec4d222f8d29dafdc22435a745..bdd3064b75167a2aa81f2f2745d0a74f165db5ba 100644
--- a/Mesh/highOrderTools.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomElastic.cpp
@@ -7,6 +7,8 @@
// Koen Hillewaert
//
+#include "OptHomElastic.h"
+#include "GModel.h"
#include "GmshConfig.h"
#if defined(HAVE_SOLVER)
@@ -21,7 +23,6 @@
#include "MPoint.h"
#include "HighOrder.h"
#include "meshGFaceOptimize.h"
-#include "highOrderTools.h"
#include "GFace.h"
#include "GRegion.h"
#include "GeomMeshMatcher.h"
@@ -36,6 +37,44 @@
#define SQU(a) ((a)*(a))
+void ElasticAnalogy(GModel *m, double threshold, bool onlyVisible)
+{
+ bool CAD, complete;
+ int meshOrder;
+
+ getMeshInfoForHighOrder(m, meshOrder, complete, CAD);
+ highOrderTools hot(m);
+ // now we smooth mesh the internal vertices of the faces
+ // we do that model face by model face
+ std::vector<MElement*> bad;
+ double worst;
+ checkHighOrderTriangles("Surface mesh", m, bad, worst);
+ {
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
+ if (onlyVisible && !(*it)->getVisibility()) continue;
+ std::vector<MElement*> v;
+ v.insert(v.begin(), (*it)->triangles.begin(), (*it)->triangles.end());
+ v.insert(v.end(), (*it)->quadrangles.begin(), (*it)->quadrangles.end());
+ if (CAD) hot.applySmoothingTo(v, (*it));
+ else hot.applySmoothingTo(v, 1.e32, false);
+ }
+ }
+ checkHighOrderTriangles("Final surface mesh", m, bad, worst);
+
+ checkHighOrderTetrahedron("Volume Mesh", m, bad, worst);
+ {
+ for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
+ if (onlyVisible && !(*it)->getVisibility())continue;
+ std::vector<MElement*> v;
+ v.insert(v.begin(), (*it)->tetrahedra.begin(), (*it)->tetrahedra.end());
+ v.insert(v.end(), (*it)->hexahedra.begin(), (*it)->hexahedra.end());
+ v.insert(v.end(), (*it)->prisms.begin(), (*it)->prisms.end());
+ hot.applySmoothingTo(v,1.e32,false);
+ }
+ }
+ checkHighOrderTetrahedron("File volume Mesh", m, bad, worst);
+}
+
void highOrderTools::moveToStraightSidedLocation(MElement *e) const
{
for(int i = 0; i < e->getNumVertices(); i++){
@@ -718,11 +757,66 @@ double highOrderTools::applySmoothingTo(std::vector<MElement*> &all,
return percentage;
}
-extern void printJacobians(GModel *m, const char *nm);
+void printJacobians(GModel *m, const char *nm)
+{
+ const int n = 100;
+ double D[n][n], X[n][n], Y[n][n], Z[n][n];
+
+ FILE *f = Fopen(nm,"w");
+ fprintf(f,"View \"\"{\n");
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ for(unsigned int j = 0; j < (*it)->triangles.size(); j++){
+ MTriangle *t = (*it)->triangles[j];
+ for(int i = 0; i < n; i++){
+ for(int k = 0; k < n - i; k++){
+ SPoint3 pt;
+ double u = (double)i / (n - 1);
+ double v = (double)k / (n - 1);
+ t->pnt(u, v, 0, pt);
+ D[i][k] = 0.; //mesh_functional_distorsion_2D(t, u, v);
+ //X[i][k] = u;
+ //Y[i][k] = v;
+ //Z[i][k] = 0.0;
+ X[i][k] = pt.x();
+ Y[i][k] = pt.y();
+ Z[i][k] = pt.z();
+ }
+ }
+ for(int i= 0; i < n -1; i++){
+ for(int k = 0; k < n - i -1; k++){
+ fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
+ X[i][k],Y[i][k],Z[i][k],
+ X[i+1][k],Y[i+1][k],Z[i+1][k],
+ X[i][k+1],Y[i][k+1],Z[i][k+1],
+ D[i][k],
+ D[i+1][k],
+ D[i][k+1]);
+ if (i != n-2 && k != n - i -2)
+ fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
+ X[i+1][k],Y[i+1][k],Z[i+1][k],
+ X[i+1][k+1],Y[i+1][k+1],Z[i+1][k+1],
+ X[i][k+1],Y[i][k+1],Z[i][k+1],
+ D[i+1][k],
+ D[i+1][k+1],
+ D[i][k+1]);
+ }
+ }
+ }
+ }
+ fprintf(f,"};\n");
+ fclose(f);
+}
void highOrderTools::makePosViewWithJacobians (const char *fn)
{
printJacobians(_gm,fn);
}
+#else
+
+void ElasticAnalogy(GModel *m, double threshold, bool onlyVisible)
+{
+ Msg::Error("Elastic analogy high-order optimzer requires the solver module");
+}
+
#endif
diff --git a/Mesh/highOrderTools.h b/contrib/HighOrderMeshOptimizer/OptHomElastic.h
similarity index 81%
rename from Mesh/highOrderTools.h
rename to contrib/HighOrderMeshOptimizer/OptHomElastic.h
index a8b7ff208eb9cebb4ad0310705945f6a51413002..15e0b83b544424b7a04b14e889da53e531165d75 100644
--- a/Mesh/highOrderTools.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomElastic.h
@@ -3,8 +3,8 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
-#ifndef _HIGH_ORDER_TOOLS_H_
-#define _HIGH_ORDER_TOOLS_H_
+#ifndef _OPT_HOM_ELASTIC_H_
+#define _OPT_HOM_ELASTIC_H_
#include <map>
#include <vector>
@@ -12,6 +12,8 @@
#include "GmshMessage.h"
#include "GModel.h"
+void ElasticAnalogy(GModel *m, double threshold, bool onlyVisible);
+
#if defined(HAVE_SOLVER)
#include "SVector3.h"
@@ -82,22 +84,6 @@ class highOrderTools
void makePosViewWithJacobians(const char *nm);
};
-#else
-
-class highOrderTools
-{
- public:
- highOrderTools(GModel *gm)
- {
- Msg::Error("Gmsh has to be compiled with solver support to use highOrderSmoother");
- }
- void applyGlobalSmoothing (){}
- void ensureMinimumDistorsion (double threshold){}
- double applySmoothingTo (GFace *gf, double tres = 0.1, bool mixed = false){ return 0.; }
- void applySmoothingTo(std::vector<MElement*> & all, GFace *gf){}
- double applySmoothingTo (std::vector<MElement*> &all, double threshold, bool mixed){ return 0.; }
-};
-
#endif
#endif
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index e8f526c1658b04793cf9376732970c4bad5741e8..4496f133e0910701838d7443a4461befbb41c220 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -9,14 +9,9 @@
#include "ParamCoord.h"
#include "polynomialBasis.h"
-
-
-
class Mesh
{
-
public:
-
Mesh(const std::set<MElement*> &els, std::set<MVertex*> & toFix, bool fixBndNodes);
inline const int &nPC() { return _nPC; }
@@ -31,13 +26,14 @@ public:
inline const int &nBezEl(int iEl) { return _nBezEl[iEl]; }
void metricMinAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
- void scaledJacAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
+ void scaledJacAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
inline int indGSJ(int iEl, int l, int iPC) { return iPC*_nBezEl[iEl]+l; }
inline double distSq(int iFV);
inline void gradDistSq(int iV, std::vector<double> &gDSq);
- void pcScale(int iFV, std::vector<double> &scale); // Calc. scale of parametric coordinates for preconditioning
+ // Calc. scale of parametric coordinates for preconditioning
+ void pcScale(int iFV, std::vector<double> &scale);
void getUvw(double *it);
void updateMesh(const double *it);
@@ -46,59 +42,69 @@ public:
void updateGEntityPositions();
void writeMSH(const char *filename);
-// inline double xyzDBG(int iV, int iCoord) { return _xyz[iV][iCoord]; }
-// inline double ixyzDBG(int iV, int iCoord) { return _ixyz[iV][iCoord]; }
-// inline double sxyzDBG(int iV, int iCoord) { return _sxyz[iV][iCoord]; }
-// inline double uvwDBG(int iFV, int iCoord) { return _uvw[iFV][iCoord]; }
-// inline int fv2VDBG(int iFV) { return _fv2V[iFV]; }
-// inline bool forcedDBG(int iV) { return _forced[iV]; }
-
private:
-
int _dim;
- int _nPC; // Total nb. of parametric coordinates
-
- std::vector<MElement*> _el; // List of elements
- std::vector<fullMatrix<double> > _scaledNormEl; // Normals to 2D elements for Jacobian regularization and scaling
- std::vector<double> _invStraightJac; // Initial Jacobians for 3D elements
- std::vector<MVertex*> _vert, _freeVert; // List of vert., free vert.
- std::vector<int> _fv2V; // Index of free vert. -> index of vert.
- std::vector<bool> _forced; // Is vertex forced?
- std::vector<SPoint3> _xyz, _ixyz; // Physical coord. of ALL vertices (current, straight, init.)
- std::vector<SPoint3> _uvw, _iuvw; // Parametric coord. of FREE vertices (current, straight, init.)
- std::vector<int> _startPCFV; // Start index of parametric coordinates for a free vertex
- std::vector<int> _nPCFV; // Number of parametric coordinates for a free vertex
- std::vector<std::vector<int> > _el2FV, _el2V; // Free vertices, vertices in element
- std::vector<int> _nBezEl, _nNodEl; // Number of Bezier poly. and nodes for an el.
- std::vector<std::vector<int> > _indPCEl; // Index of parametric coord. for an el.
+ // Total nb. of parametric coordinates
+ int _nPC;
+ // List of elements
+ std::vector<MElement*> _el;
+ // Normals to 2D elements for Jacobian regularization and scaling
+ std::vector<fullMatrix<double> > _scaledNormEl;
+ // Initial Jacobians for 3D elements
+ std::vector<double> _invStraightJac;
+ // List of vert., free vert.
+ std::vector<MVertex*> _vert, _freeVert;
+ // Index of free vert. -> index of vert.
+ std::vector<int> _fv2V;
+ // Is vertex forced?
+ std::vector<bool> _forced;
+ // Physical coord. of ALL vertices (current, straight, init.)
+ std::vector<SPoint3> _xyz, _ixyz;
+ // Parametric coord. of FREE vertices (current, straight, init.)
+ std::vector<SPoint3> _uvw, _iuvw;
+ // Start index of parametric coordinates for a free vertex
+ std::vector<int> _startPCFV;
+ // Number of parametric coordinates for a free vertex
+ std::vector<int> _nPCFV;
+ // Free vertices, vertices in element
+ std::vector<std::vector<int> > _el2FV, _el2V;
+ // Number of Bezier poly. and nodes for an el.
+ std::vector<int> _nBezEl, _nNodEl;
+ // Index of parametric coord. for an el.
+ std::vector<std::vector<int> > _indPCEl;
ParamCoord *_pc;
int addVert(MVertex* vert);
- int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);
+ int addFreeVert(MVertex* vert, const int iV, const int nPCV,
+ std::set<MVertex*> &toFix);
void calcScaledNormalEl2D(int iEl);
- static inline int indJB2DBase(int nNod, int l, int i, int j) { return (l*nNod+i)*nNod+j; }
- inline int indJB2D(int iEl, int l, int i, int j) { return indJB2DBase(_nNodEl[iEl],l,i,j); }
- static inline int indJB3DBase(int nNod, int l, int i, int j, int m) { return ((l*nNod+i)*nNod+j)*nNod+m; }
- inline int indJB3D(int iEl, int l, int i, int j, int m) { return indJB3DBase(_nNodEl[iEl],l,i,j,m); }
-
+ static inline int indJB2DBase(int nNod, int l, int i, int j)
+ {
+ return (l*nNod+i)*nNod+j;
+ }
+ inline int indJB2D(int iEl, int l, int i, int j)
+ {
+ return indJB2DBase(_nNodEl[iEl],l,i,j);
+ }
+ static inline int indJB3DBase(int nNod, int l, int i, int j, int m)
+ {
+ return ((l*nNod+i)*nNod+j)*nNod+m;
+ }
+ inline int indJB3D(int iEl, int l, int i, int j, int m)
+ {
+ return indJB3DBase(_nNodEl[iEl],l,i,j,m);
+ }
};
-
-
double Mesh::distSq(int iFV)
{
-
SPoint3 d = _xyz[_fv2V[iFV]]-_ixyz[_fv2V[iFV]];
return d[0]*d[0]+d[1]*d[1]+d[2]*d[2];
-
}
-
-
void Mesh::gradDistSq(int iFV, std::vector<double> &gDSq)
{
-
SPoint3 gXyz = (_xyz[_fv2V[iFV]]-_ixyz[_fv2V[iFV]])*2.;
SPoint3 gUvw;
_pc->gXyz2gUvw(_freeVert[iFV],_uvw[iFV],gXyz,gUvw);
@@ -106,9 +112,6 @@ void Mesh::gradDistSq(int iFV, std::vector<double> &gDSq)
gDSq[0] = gUvw[0];
if (_nPCFV[iFV] >= 2) gDSq[1] = gUvw[1];
if (_nPCFV[iFV] == 3) gDSq[2] = gUvw[2];
-
}
-
-
#endif
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 7ee49c94b42f9240acc9cc120d8d5f864d385233..940dfd3b4c98f9e18ce6d84cdb12c8af48377273 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -11,17 +11,11 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MLine.h"
-#include "highOrderTools.h"
#include "OptHomMesh.h"
#include "OptHOM.h"
#include "OS.h"
#include <stack>
-#ifdef HAVE_FLTK
-#include "highOrderToolsWindow.h"
-#include "FlGui.h"
-#endif
-
double distMaxStraight(MElement *el)
{
const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();