diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index d715c6f93ebfdde34d5032b0d524f93c16248abf..de16bb1b7e08e41fc7cf587e1ec499286c247d49 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -452,7 +452,7 @@ double GFace::curvatureDiv(const SPoint2 ¶m) const
const double eps = 1.e-5;
- Pair<SVector3,SVector3> der = firstDer(param);
+ Pair<SVector3, SVector3> der = firstDer(param);
SVector3 du = der.first();
SVector3 dv = der.second();
@@ -463,18 +463,20 @@ double GFace::curvatureDiv(const SPoint2 ¶m) const
du.normalize();
dv.normalize();
- SVector3 n1,n2,n3,n4;
- if( param.x() - eps < 0.0 ) {
+ SVector3 n1, n2, n3, n4;
+ if(param.x() - eps < 0.0) {
n1 = normal(SPoint2(param.x(), param.y()));
n2 = normal(SPoint2(param.x() + eps, param.y()));
- } else {
+ }
+ else {
n1 = normal(SPoint2(param.x() - eps, param.y()));
n2 = normal(SPoint2(param.x(), param.y()));
}
- if( param.y() - eps < 0.0 ) {
+ if(param.y() - eps < 0.0) {
n3 = normal(SPoint2(param.x(), param.y() ));
n4 = normal(SPoint2(param.x(), param.y() + eps));
- } else {
+ }
+ else {
n3 = normal(SPoint2(param.x(), param.y() - eps));
n4 = normal(SPoint2(param.x(), param.y() ));
}
@@ -485,12 +487,12 @@ double GFace::curvatureDiv(const SPoint2 ¶m) const
SVector3 dudu = SVector3();
SVector3 dvdv = SVector3();
SVector3 dudv = SVector3();
- secondDer(param,&dudu,&dvdv,&dudv);
+ secondDer(param, &dudu, &dvdv, &dudv);
- double ddu = dot(dndu,du);
- double ddv = dot(dndv,dv);
+ double ddu = dot(dndu, du);
+ double ddv = dot(dndv, dv);
- return ( fabs(ddu) + fabs(ddv) ) / detJ;
+ return (fabs(ddu) + fabs(ddv)) / detJ;
}
double GFace::curvatureMax(const SPoint2 ¶m) const
@@ -503,31 +505,26 @@ double GFace::curvatureMax(const SPoint2 ¶m) const
return fabs(eigVal[1]);
}
-double GFace::curvatures(const SPoint2 ¶m,
- SVector3 *dirMax,
- SVector3 *dirMin,
- double *curvMax,
- double *curvMin) const
+double GFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const
{
- Pair<SVector3,SVector3> D1 = firstDer(param);
+ Pair<SVector3, SVector3> D1 = firstDer(param);
- if (geomType() == Plane)
- {
- *dirMax = D1.first();
- *dirMin = D1.second();
- *curvMax = 0.;
- *curvMin = 0.;
- return 0.;
- }
-
- if (geomType() == Sphere)
- {
- *dirMax = D1.first();
- *dirMin = D1.second();
- *curvMax = curvatureDiv(param);
- *curvMin = *curvMax;
- return *curvMax;
- }
+ if(geomType() == Plane){
+ *dirMax = D1.first();
+ *dirMin = D1.second();
+ *curvMax = 0.;
+ *curvMin = 0.;
+ return 0.;
+ }
+
+ if(geomType() == Sphere){
+ *dirMax = D1.first();
+ *dirMin = D1.second();
+ *curvMax = curvatureDiv(param);
+ *curvMin = *curvMax;
+ return *curvMax;
+ }
double eigVal[2], eigVec[8];
getMetricEigenVectors(param, eigVal, eigVec);
@@ -550,8 +547,7 @@ double GFace::getMetricEigenvalue(const SPoint2 &)
// eigen values are absolute values and sorted from min to max of absolute values
// eigen vectors are the corresponding COLUMNS of eigVec
void GFace::getMetricEigenVectors(const SPoint2 ¶m,
- double eigVal[2],
- double eigVec[4]) const
+ double eigVal[2], double eigVec[4]) const
{
// first derivatives
Pair<SVector3,SVector3> D1 = firstDer(param);
@@ -563,7 +559,7 @@ void GFace::getMetricEigenVectors(const SPoint2 ¶m,
SVector3 dudu = SVector3();
SVector3 dvdv = SVector3();
SVector3 dudv = SVector3();
- secondDer(param,&dudu,&dvdv,&dudv);
+ secondDer(param, &dudu, &dvdv, &dudv);
// first form
double form1[2][2];
@@ -579,7 +575,7 @@ void GFace::getMetricEigenVectors(const SPoint2 ¶m,
// inverse of first form
double inv_form1[2][2];
- double inv_det_form1 = 1. / ( form1[0][0] * form1[1][1] - form1[1][0] * form1[0][1] );
+ double inv_det_form1 = 1. / (form1[0][0] * form1[1][1] - form1[1][0] * form1[0][1]);
inv_form1[0][0] = inv_det_form1 * form1[1][1];
inv_form1[1][1] = inv_det_form1 * form1[0][0];
inv_form1[1][0] = inv_form1[0][1] = -1 * inv_det_form1 * form1[0][1];
@@ -595,17 +591,17 @@ void GFace::getMetricEigenVectors(const SPoint2 ¶m,
int work1[2];
double work2[2];
double eigValI[2];
- if ( EigSolve(2,2,N,eigVal,eigValI,eigVec,work1,work2) != 1 ) {
- Msg::Warning("Problem in eigen vectors computation");
- printf(" N: %f %f %f %f\n",N[0],N[1],N[2],N[3]);
- printf(" * Eigen values:\n %f + i * %f\n %f + i * %f\n",
- eigVal[0],eigValI[0],eigVal[1],eigValI[1]);
- printf(" * Eigen vectors (trust it only if eigen values are real):\n");
- printf(" ( %f, %f ),\n ( %f, %f ).\n",
- eigVec[0],eigVec[2],eigVec[1],eigVec[3]);
- throw;
- }
- if ( fabs(eigValI[0]) > 1.e-12 || fabs(eigValI[1]) > 1.e-12 ) {
+ if (EigSolve(2, 2, N, eigVal, eigValI, eigVec, work1, work2) != 1) {
+ Msg::Error("Problem in eigen vectors computation");
+ Msg::Error(" N: %f %f %f %f", N[0], N[1], N[2], N[3]);
+ Msg::Error(" * Eigen values:");
+ Msg::Error(" %f + i * %f, %f + i * %f",
+ eigVal[0], eigValI[0], eigVal[1], eigValI[1]);
+ Msg::Error(" * Eigen vectors (trust it only if eigen values are real):");
+ Msg::Error(" ( %f, %f ), ( %f, %f )",
+ eigVec[0], eigVec[2], eigVec[1], eigVec[3]);
+ }
+ if (fabs(eigValI[0]) > 1.e-12 || fabs(eigValI[1]) > 1.e-12) {
Msg::Error("Found imaginary eigenvalues");
}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index 814bb3db4589678a6b8157577e7fe6ad7307e3a8..ad40f404869911b81c064cbca379f5fa4f82c6d6 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -148,10 +148,9 @@ class GFace : public GEntity
virtual SVector3 normal(const SPoint2 ¶m) const;
// return the first derivate of the face at the parameter location
- virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const = 0;
+ virtual Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const = 0;
// compute the second derivates of the face at the parameter location
- // (default implementation by central differences)
// the derivates have to be allocated before calling this function
virtual void secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const = 0;
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index bfb3493b3e6037332cb46e96bde1a6f4011f247e..fdf1e6c1489a823a3f3ffd0ac7e9825230f20a31 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -141,7 +141,7 @@ GPoint OCCFace::point(double par1, double par2) const
return GPoint(val.X(), val.Y(), val.Z(), this, pp);
}
-GPoint OCCFace::closestPoint(const SPoint3 & qp, const double initialGuess[2]) const
+GPoint OCCFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) const
{
gp_Pnt pnt(qp.x(), qp.y(), qp.z());
GeomAPI_ProjectPointOnSurf proj(pnt, occface, umin, umax, vmin, vmax);
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 3aadb7f106c7996918fcfe798447c087cbd015ab..bf40406911937f28a172d0ffad54c9ab8b81b40f 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -26,13 +26,10 @@ discreteEdge::discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1)
Tree_Add(model->getGEOInternals()->Curves, &c);
CreateReversedCurve(c);
#endif
-
}
-
void discreteEdge::orderMLines()
{
-
//printf(" *** ORDERING DISCRETE EDGE %d of size %d \n", this->tag(), lines.size());
std::vector<MLine*> _m ;
@@ -281,11 +278,9 @@ void discreteEdge::parametrize()
}
-void discreteEdge::getLocalParameter ( const double &t,
- int &iLine,
- double & tLoc) const
+void discreteEdge::getLocalParameter(const double &t, int &iLine,
+ double &tLoc) const
{
-
for (iLine=0 ; iLine<lines.size() ;iLine++){
double tmin = _pars[iLine];
double tmax = _pars[iLine+1];
@@ -319,7 +314,6 @@ GPoint discreteEdge::point(double par) const
SVector3 discreteEdge::firstDer(double par) const
{
-
double tLoc;
int iEdge;
getLocalParameter(par,iEdge,tLoc);
@@ -335,11 +329,9 @@ SVector3 discreteEdge::firstDer(double par) const
//printf("firstDer discreteEdge par=%g, dx=%g, dy=%g dz=%g dt=%g \n", par,dx,dy,dz, dt);
return SVector3(dx,dy,dz);
-
}
-Range<double> discreteEdge::parBounds(int i) const {
-
- return Range<double>(0, lines.size());
-
+Range<double> discreteEdge::parBounds(int i) const
+{
+ return Range<double>(0, lines.size());
}
diff --git a/Geo/discreteEdge.h b/Geo/discreteEdge.h
index a0cc4eac2b131c520b94200bb6935469b8998d8f..1e6d131a249c99dc2edc8188ecdce8555e8424d0 100644
--- a/Geo/discreteEdge.h
+++ b/Geo/discreteEdge.h
@@ -18,15 +18,13 @@ class discreteEdge : public GEdge {
public:
discreteEdge(GModel *model, int num, GVertex *_v0, GVertex *_v1);
virtual ~discreteEdge() {}
- void getLocalParameter ( const double &t,
- int &iEdge,
- double & tLoc) const;
+ void getLocalParameter(const double &t, int &iEdge, double &tLoc) const;
virtual GeomType geomType() const { return DiscreteCurve; }
virtual GPoint point(double p) const;
virtual SVector3 firstDer(double par) const;
virtual Range<double> parBounds(int) const;
- void parametrize() ;
- void orderMLines() ;
+ void parametrize();
+ void orderMLines();
void setBoundVertices();
};
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index 8de241977ef7547f59ed032c9b414a962a1f90df..439a7ec1bb18b2df8bb3c59b556db882fa8fd051 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -22,8 +22,8 @@ discreteFace::discreteFace(GModel *model, int num) : GFace(model, num)
meshStatistics.status = GFace::DONE;
}
-void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges){
-
+void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges)
+{
//find the boundary edges
std::list<MEdge> bound_edges;
for (int iFace = 0; iFace < getNumMeshElements() ; iFace++) {
@@ -39,7 +39,6 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_
}
}
-
//for the boundary edges, associate the tag of the current discrete face
for (std::list<MEdge>::iterator itv = bound_edges.begin() ; itv != bound_edges.end() ; ++itv){
std::map<MEdge, std::vector<int> , Less_Edge >::iterator itmap = map_edges.find(*itv);
@@ -57,8 +56,6 @@ void discreteFace::findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_
}
//printf( "There are %d bound msh edges \n ", map_edges.size());
-
-
}
void discreteFace::setBoundEdges(std::vector<int> tagEdges)
@@ -85,11 +82,8 @@ void discreteFace::setBoundEdges(std::vector<int> tagEdges)
}
// printf("bound edges =%d \n", edges().size());
-
-
}
-
GPoint discreteFace::point(double par1, double par2) const
{
Msg::Error("Cannot evaluate point on discrete face");
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index d44f1ca8fc385caedc23644192b469b5940f2b51..87ed1db371b2153eefba5257e1fe3318e3e3baab 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -19,11 +19,11 @@ class discreteFace : public GFace {
virtual SPoint2 parFromPoint(const SPoint3 &p) const;
virtual SVector3 normal(const SPoint2 ¶m) const;
virtual GEntity::GeomType geomType() const { return DiscreteSurface; }
- virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ virtual Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const;
- void setBoundEdges( std::vector<int> tagEdges );
- void findEdges( std::map<MEdge, std::vector<int>, Less_Edge > &map_edges);
+ void setBoundEdges(std::vector<int> tagEdges);
+ void findEdges(std::map<MEdge, std::vector<int>, Less_Edge > &map_edges);
};
#endif
diff --git a/Geo/fourierFace.cpp b/Geo/fourierFace.cpp
index 4203b1acfff3ae25441845a2a17a6cf1cff4657c..da911c94f4bb376790fd9e51fae9875484f0c01e 100644
--- a/Geo/fourierFace.cpp
+++ b/Geo/fourierFace.cpp
@@ -74,7 +74,7 @@ Pair<SVector3, SVector3> fourierFace::firstDer(const SPoint2 ¶m) const
void fourierFace::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
{
- Msg::Error("Computation of the second derivatives not implemented for fourer face");
+ Msg::Error("Computation of the second derivatives not implemented for fourier face");
}
#endif
diff --git a/Geo/fourierProjectionFace.cpp b/Geo/fourierProjectionFace.cpp
index a8a9001c115ecd17a071fbe86133e2854afe2e2d..3d9174209e3aa64369e21f3afb741e5a90cb45c6 100644
--- a/Geo/fourierProjectionFace.cpp
+++ b/Geo/fourierProjectionFace.cpp
@@ -4,6 +4,7 @@
// bugs and problems to <gmsh@geuz.org>.
#include "GmshConfig.h"
+#include "GmshMessage.h"
#include "fourierProjectionFace.h"
#include "VertexArray.h"
@@ -40,6 +41,12 @@ Pair<SVector3,SVector3> fourierProjectionFace::firstDer(const SPoint2 ¶m) co
return Pair<SVector3,SVector3>(du,dv);
}
+void fourierProjectionFace::secondDer(const SPoint2 ¶m,
+ SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
+{
+ Msg::Error("Computation of the second derivatives not implemented");
+}
+
SVector3 fourierProjectionFace::normal(const SPoint2 ¶m) const
{
double x, y, z;
diff --git a/Geo/fourierProjectionFace.h b/Geo/fourierProjectionFace.h
index e27a1713ff732eaaa2f3068ec5f7cc5d9e979519..95adf83a65873fdec6bf1f2aae400b17fddd6237 100644
--- a/Geo/fourierProjectionFace.h
+++ b/Geo/fourierProjectionFace.h
@@ -24,6 +24,7 @@ class fourierProjectionFace : public GFace {
GPoint point(double par1, double par2) const;
SVector3 normal(const SPoint2 ¶m) const;
Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
+ void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
SPoint2 parFromPoint(const SPoint3 &) const;
virtual GEntity::GeomType geomType() const { return GEntity::ProjectionFace; }
ModelType getNativeType() const { return UnknownModel; }