diff --git a/CMakeLists.txt b/CMakeLists.txt
index 44a647a61ce0b9dd112db0f6c2207f5219b6a116..51e82388f165878e668d9bc5d6fbf810d4e636aa 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -18,7 +18,7 @@ project(gmsh CXX C)
option(ENABLE_ACIS "Enable ACIS geometrical models" ON)
option(ENABLE_ANN "Enable ANN to compute Approximate Nearest Neighbors" ON)
-option(ENABLE_APP_BUNDLE_PATH "Use special install path for .app bundle on Mac" ON)
+option(ENABLE_APP_BUNDLE "Create .app bundle on Mac when installing" ON)
option(ENABLE_BAMG "Enable Bamg mesh generator" ON)
option(ENABLE_BFGS "Enable BFGS" ON)
option(ENABLE_BLAS_LAPACK "Use BLAS and Lapack for linear algebra" ON)
@@ -48,7 +48,7 @@ option(ENABLE_PARSER "Build the GEO file parser" ON)
option(ENABLE_PETSC "Enable PETSc linear algebra solvers" ON)
option(ENABLE_PLUGINS "Build the post-processing plugins" ON)
option(ENABLE_POST "Build the post-processing module" ON)
-option(ENABLE_QT "Build QT GUI (not functional: for testing only)" OFF)
+option(ENABLE_QT "Build proof-of-concept QT GUI" OFF)
option(ENABLE_RBF "Enable RBF project" OFF)
option(ENABLE_SALOME "Enable Salome routines for CAD healing" ON)
option(ENABLE_SLEPC "Enable SLEPc eigensolvers" ON)
@@ -135,8 +135,7 @@ include(CheckFunctionExists)
include(CheckIncludeFile)
if(MSVC)
- # remove really annoying (and stupid, and wrong) warning about
- # bool/int cast performance in Visual C++
+ # remove annoying warning about bool/int cast performance
set(GMSH_CONFIG_PRAGMAS "#pragma warning(disable:4800)")
if(ENABLE_MSVC_STATIC_RUNTIME)
foreach(VAR
@@ -1046,7 +1045,7 @@ if(WIN32 OR CYGWIN)
unix2dos(TUTORIAL_FILES)
unix2dos(DEMO_FILES)
endif(CYGWIN)
-elseif(APPLE AND ENABLE_APP_BUNDLE_PATH)
+elseif(APPLE AND ENABLE_APP_BUNDLE)
# set these so that the files get installed nicely in the MacOSX
# .app bundle
set(GMSH_BIN ../MacOS)
@@ -1170,7 +1169,7 @@ if(UNIX)
set(CPACK_INSTALL_COMMANDS "rm -rf ${CMAKE_CURRENT_BINARY_DIR}/_CPack_Packages")
endif(UNIX)
-if(APPLE)
+if(APPLE AND ENABLE_APP_BUNDLE)
set(CPACK_GENERATOR Bundle)
set(CPACK_BUNDLE_NAME Gmsh)
file(READ ${CMAKE_CURRENT_SOURCE_DIR}/utils/misc/gmsh_app.plist F0)
@@ -1178,20 +1177,18 @@ if(APPLE)
file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/Info.plist "${F1}")
set(CPACK_BUNDLE_PLIST ${CMAKE_CURRENT_BINARY_DIR}/Info.plist)
set(CPACK_BUNDLE_ICON ${CMAKE_CURRENT_SOURCE_DIR}/Fltk/MacIcons.icns)
- if(ENABLE_APP_BUNDLE_PATH)
- install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/Fltk/MacIconsGeo.icns DESTINATION .
- RENAME GmshGeo.icns)
- install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/Fltk/MacIconsMsh.icns DESTINATION .
- RENAME GmshMsh.icns)
- install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/Fltk/MacIconsPos.icns DESTINATION .
- RENAME GmshPos.icns)
- endif(ENABLE_APP_BUNDLE_PATH)
+ install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/Fltk/MacIconsGeo.icns DESTINATION .
+ RENAME GmshGeo.icns)
+ install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/Fltk/MacIconsMsh.icns DESTINATION .
+ RENAME GmshMsh.icns)
+ install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/Fltk/MacIconsPos.icns DESTINATION .
+ RENAME GmshPos.icns)
set(CPACK_PACKAGE_ICON ${CMAKE_CURRENT_SOURCE_DIR}/Fltk/MacIcons.icns)
elseif(WIN32 OR CYGWIN)
set(CPACK_GENERATOR ZIP)
-else(APPLE)
+else(APPLE AND ENABLE_APP_BUNDLE)
set(CPACK_GENERATOR TGZ)
-endif(APPLE)
+endif(APPLE AND ENABLE_APP_BUNDLE)
include(CPack)
@@ -1214,4 +1211,3 @@ message("")
mark_as_advanced(BISON FLEX GMP_LIB GMSH_EXTRA_VERSION HDF5_LIB MAKEINFO
MED_LIB OCC_INC SZ_LIB TAUCS_LIB ACIS_LIB TEXI2PDF)
-
diff --git a/Common/Context.h b/Common/Context.h
index f76ebcb882c73ed48b13fbba8b5cea8021f8900f..ab570ec0080f58d69f0b4ef8930069f646500dc4 100644
--- a/Common/Context.h
+++ b/Common/Context.h
@@ -55,8 +55,6 @@ struct contextGeometryOptions {
int matchGeomAndMesh;
};
-
-
class CTX {
private:
static CTX *_instance;
diff --git a/Common/OctreeInternals.cpp b/Common/OctreeInternals.cpp
index b2b4fc8a3be59a8f4d186df4e03bd6a807d20c60..bfd2d2fca8345c063da49ff58d0c8fa5bdee2c26 100644
--- a/Common/OctreeInternals.cpp
+++ b/Common/OctreeInternals.cpp
@@ -221,7 +221,7 @@ octantBucket *findElementBucket(octantBucket *_buckets_head, double *_pt)
}
} // for loop i
if (i == num) {
- // printf("Error, no bucket contains the given point! ");
+ // Msg::Error("No bucket contains the given point!");
return NULL;
}
} // for while loop
@@ -429,6 +429,6 @@ void *searchAllElements(octantBucket *_buckets_head, double *_pt, globalInfo *_g
if (flag1)
return (void *)(_elements);
- // Msg::Warning("This point is not found in any element! It is not in the domain");
+ // Msg::Warning("This point is not found in any element! It is not in the domain");
return NULL;
}
diff --git a/Geo/ACISEdge.cpp b/Geo/ACISEdge.cpp
index f9bf010439d8b545edceaa0a1d4e71914a4b2e14..c9444b51621e3b44b6ed4e61fe0e5607f5599369 100644
--- a/Geo/ACISEdge.cpp
+++ b/Geo/ACISEdge.cpp
@@ -35,7 +35,6 @@ GEdge *getACISEdgeByNativePtr(GModel *model, EDGE * toFind)
return 0;
}
-
ACISEdge::ACISEdge(GModel *model, EDGE* edge, int num, GVertex *v1, GVertex *v2)
: GEdge(model, num, v1, v2), _e(edge)
{
@@ -49,8 +48,6 @@ ACISEdge::ACISEdge(GModel *model, EDGE* edge, int num, GVertex *v1, GVertex *v2)
// s1 += _cur->equation().param_period()/2;
}
Msg::Info("ACIS Edge %d is %s goes from %g to %g",tag(),getTypeString().c_str(),s0,s1);
-
-
}
Range<double> ACISEdge::parBounds(int i) const
@@ -68,11 +65,11 @@ SPoint2 ACISEdge::reparamOnFace(const GFace *face, double epar, int dir) const
return pt2;
}
-GPoint ACISEdge::closestPoint(const SPoint3 &qp, double ¶m) const{
+GPoint ACISEdge::closestPoint(const SPoint3 &qp, double ¶m) const
+{
throw;
}
-
// True if the edge is a seam for the given face
bool ACISEdge::isSeam(const GFace *face) const
{
@@ -84,7 +81,8 @@ bool ACISEdge::isSeam(const GFace *face) const
!(((FACE*)face->getNativePtr())->geometry()->equation().periodic_v()))
{
SPAinterval cur_rang = _cur->equation().param_range();
- SPAinterval sur_rang_u = ((FACE*)face->getNativePtr())->geometry()->equation().param_range_u();
+ SPAinterval sur_rang_u =
+ ((FACE*)face->getNativePtr())->geometry()->equation().param_range_u();
SPAposition p1, p2;
p1 = _cur->equation().eval_position(cur_rang.start_pt());
p2 = _cur->equation().eval_position(cur_rang.end_pt());
@@ -102,7 +100,8 @@ bool ACISEdge::isSeam(const GFace *face) const
(((FACE*)face->getNativePtr())->geometry()->equation().periodic_v()))
{
SPAinterval cur_rang = _cur->equation().param_range();
- SPAinterval sur_rang_v = ((FACE*)face->getNativePtr())->geometry()->equation().param_range_v();
+ SPAinterval sur_rang_v =
+ ((FACE*)face->getNativePtr())->geometry()->equation().param_range_v();
SPAposition p1, p2;
p1 = _cur->equation().eval_position(cur_rang.start_pt());
p2 = _cur->equation().eval_position(cur_rang.end_pt());
@@ -120,8 +119,10 @@ bool ACISEdge::isSeam(const GFace *face) const
(((FACE*)face->getNativePtr())->geometry()->equation().periodic_v()))
{
SPAinterval cur_rang = _cur->equation().param_range();
- SPAinterval sur_rang_u = ((FACE*)face->getNativePtr())->geometry()->equation().param_range_u();
- SPAinterval sur_rang_v = ((FACE*)face->getNativePtr())->geometry()->equation().param_range_v();
+ SPAinterval sur_rang_u =
+ ((FACE*)face->getNativePtr())->geometry()->equation().param_range_u();
+ SPAinterval sur_rang_v =
+ ((FACE*)face->getNativePtr())->geometry()->equation().param_range_v();
SPAposition p1, p2;
p1 = _cur->equation().eval_position(cur_rang.start_pt());
p2 = _cur->equation().eval_position(cur_rang.end_pt());
@@ -139,11 +140,11 @@ bool ACISEdge::isSeam(const GFace *face) const
}
}
-bool ACISEdge::degenerate(int) const{
+bool ACISEdge::degenerate(int) const
+{
return _e->length() == 0.0;
}
-
GPoint ACISEdge::point(double par) const
{
CURVE *c = _e->geometry();
diff --git a/Geo/ACISFace.cpp b/Geo/ACISFace.cpp
index 09d4b6ef817c8592c78e7c9ca414b29b4d560625..0eb6738e17679b55fea0d359e29ad984e7a43f5e 100644
--- a/Geo/ACISFace.cpp
+++ b/Geo/ACISFace.cpp
@@ -71,12 +71,10 @@ ACISFace::ACISFace(GModel *m, FACE *f, int num)
vmin = _v.start_pt();
vmax = _v.end_pt();
-
if ( surf.closed_u() || surf.closed_v() ) {
printf("%d %g %g -- %g %g\n",tag(),umin,umax,vmin,vmax);
}
-
_periodic[0] = surf.periodic_u();
_periodic[1] = surf.periodic_v();
@@ -97,7 +95,6 @@ Range<double> ACISFace::parBounds(int i) const
SVector3 ACISFace::normal(const SPoint2 ¶m) const
{
-
SURFACE *Surf = _f->geometry();
const surface &surf = Surf->equation();
@@ -112,10 +109,11 @@ Pair<SVector3,SVector3> ACISFace::firstDer(const SPoint2 ¶m) const
SPAposition pos;
SPAvector first_derivs[2];
surf.eval (SPApar_pos(param.x(),param.y()), pos,first_derivs);
- return Pair<SVector3,SVector3> (
- SVector3(first_derivs[0].component(0),first_derivs[0].component(1),first_derivs[0].component(2)),
- SVector3(first_derivs[1].component(0),first_derivs[1].component(1),first_derivs[1].component(2))
- );
+ return Pair<SVector3,SVector3>
+ (SVector3(first_derivs[0].component(0), first_derivs[0].component(1),
+ first_derivs[0].component(2)),
+ SVector3(first_derivs[1].component(0), first_derivs[1].component(1),
+ first_derivs[1].component(2)));
}
void ACISFace::secondDer(const SPoint2 ¶m,
@@ -153,12 +151,8 @@ SPoint2 ACISFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
SURFACE *Surf = _f->geometry();
SPAposition pt(qp.x(),qp.y(),qp.z());
SPApar_pos ppt = Surf->equation().param(pt);
-
- SPoint2 pt2(ppt.u,ppt.v); // 200306
-
+ SPoint2 pt2(ppt.u,ppt.v);
GPoint sp = point(ppt.u,ppt.v);
- // printf("%g %g %g vs %g %g %g\n",qp.x(),qp.y(),qp.z(),sp.x(),sp.y(),sp.z());
-
return pt2;
}
@@ -203,7 +197,6 @@ double ACISFace::period(int dir) const {
return _f->geometry()->equation().param_period_v();
}
-
bool ACISFace::buildSTLTriangulation(bool force)
{
if(stl_triangles.size()){
@@ -231,5 +224,4 @@ GFace *getACISFaceByNativePtr(GModel *model, FACE *f)
return 0;
}
-
#endif
diff --git a/Geo/CGNSOptions.h b/Geo/CGNSOptions.h
index ee6bf182029dddec88b5c99b5205695d8ccf819d..4d2823a5394ff7243523516cfb8e6b39a6b11ca6 100644
--- a/Geo/CGNSOptions.h
+++ b/Geo/CGNSOptions.h
@@ -11,18 +11,15 @@
class CGNSOptions
{
public:
- // Types
enum CGNSLocationType {
LocVertex = 0,
LocFace = 1
};
- // Data
std::string baseName;
std::string zoneName;
std::string interfaceName;
std::string patchName;
-
int gridConnectivityLocation; // Location of connectivity(values
// CGNSLocationType)
int bocoLocation; // Location of BC (values
@@ -37,15 +34,10 @@ class CGNSOptions
bool writeUserDef; // T - write user-defined elements for
// element types unsupported by CGNS
-//--Constructor
-
CGNSOptions()
{
setDefaults();
}
-
-//--Default values
-
void setDefaults()
{
baseName = "Base_1";
diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index c2c147c0125473fe41e516aff98279cd0d82415c..9480e36f19edc4bf7b027f0f7ca45f5d309595f3 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -1,7 +1,13 @@
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by Matti Pellikka <matti.pellikka@tut.fi>.
+
#include "CellComplex.h"
#include "MElement.h"
-
CellComplex::CellComplex(GModel* model,
std::vector<MElement*>& domainElements,
std::vector<MElement*>& subdomainElements) :
@@ -11,7 +17,6 @@ CellComplex::CellComplex(GModel* model,
_insertCells(subdomainElements, 1);
_insertCells(domainElements, 0);
-
int num = 0;
for(int dim = 0; dim < 4; dim++){
if(getSize(dim) != 0) _dim = dim;
diff --git a/Geo/CellComplex.h b/Geo/CellComplex.h
index cebc1567fb7d11017fa3634917a4f2cd88722f7c..45c9dad5ba1b0af3684b768fc208ebe7a0f51781 100644
--- a/Geo/CellComplex.h
+++ b/Geo/CellComplex.h
@@ -21,10 +21,8 @@
class Cell;
class BdInfo;
-
class CellComplex
{
-
private:
GModel* _model;
@@ -136,7 +134,6 @@ class CellComplex
// experimental
int saveComplex(std::string filename);
int loadComplex(std::string filename);
-
};
#endif
diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 99f2362b444f7c7af300c48337819a4a2b3df2d8..20cae44ffb0288644e6d1e728e5ba2c07844a5b2 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -3,6 +3,9 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
+#include <iostream>
+#include <fstream>
+#include <cmath>
#include "Curvature.h"
#include "MElement.h"
#include "MTriangle.h"
@@ -13,50 +16,25 @@
#include "OS.h"
#include "SBoundingBox3d.h"
-#include<iostream>
-#include<fstream>
-#include<cmath>
-
#define NEXT(i) ((i)<2 ? (i)+1 : (i)-2)
#define PREV(i) ((i)>0 ? (i)-1 : (i)+2)
-//========================================================================================================
-//Initialization of the static variables:
Curvature* Curvature::_instance = 0;
bool Curvature::_destroyed = false;
bool Curvature::_alreadyComputedCurvature = false;
-//========================================================================================================
-
-//CONSTRUCTOR
-Curvature::Curvature(){
-}
-
-// Curvature::Curvature(std::list<GFace*> &myFaces){
-
-// std::list<GFace*>::const_iterator it = myFaces.begin();
-// for( ; it != myFaces.end() ; ++it){
-// _ptFinalEntityList.push_back(*it);
-// }
-
-// }
-
-//========================================================================================================
-
Curvature::~Curvature()
{
- _instance = 0;
+ _instance = 0;
_destroyed = true;
-
}
-//========================================================================================================
+
void Curvature::onDeadReference()
{
std::cout << "Dead reference of Curvature detected" << std::endl;
}
-//========================================================================================================
-Curvature& Curvature::getInstance()
+Curvature &Curvature::getInstance()
{
if (!_instance) {
if(_destroyed){
@@ -68,74 +46,63 @@ Curvature& Curvature::getInstance()
}
return *_instance;
}
-//========================================================================================================
- bool Curvature::valueAlreadyComputed()
- {
- return _alreadyComputedCurvature;
- }
-//========================================================================================================
+bool Curvature::valueAlreadyComputed()
+{
+ return _alreadyComputedCurvature;
+}
- void Curvature::create()
- {
- static Curvature instance;
- _instance = & instance;
- }
+void Curvature::create()
+{
+ static Curvature instance;
+ _instance = & instance;
+}
- //========================================================================================================
- void Curvature::retrieveCompounds() {
+void Curvature::retrieveCompounds()
+{
#if defined(HAVE_SOLVER)
- std::vector<GEntity*> entities;
- _model->getEntities(entities);
-
- for(int ie = 0; ie < entities.size(); ++ie) {
-
- if( entities[ie]->geomType() == GEntity::CompoundSurface ) {
- GFaceCompound* compound = dynamic_cast<GFaceCompound*>(entities[ie]);
- std::list<GFace*> tempcompounds = compound->getCompounds();
- std::list<GFace*>::iterator surfIterator;
-
- for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end(); ++surfIterator) {
- if ((*surfIterator)->geomType() == GEntity::DiscreteSurface) {
- _ptFinalEntityList.push_back(*surfIterator);
- }
- }
- }
- else if (entities[ie]->geomType() == GEntity::DiscreteSurface) {
- _ptFinalEntityList.push_back(dynamic_cast<GFace*>(entities[ie]));
- }
- }
+ std::vector<GEntity*> entities;
+ _model->getEntities(entities);
+
+ for(int ie = 0; ie < entities.size(); ++ie) {
+
+ if(entities[ie]->geomType() == GEntity::CompoundSurface) {
+ GFaceCompound* compound = dynamic_cast<GFaceCompound*>(entities[ie]);
+ std::list<GFace*> tempcompounds = compound->getCompounds();
+ std::list<GFace*>::iterator surfIterator;
+ for(surfIterator = tempcompounds.begin(); surfIterator != tempcompounds.end();
+ ++surfIterator) {
+ if ((*surfIterator)->geomType() == GEntity::DiscreteSurface) {
+ _ptFinalEntityList.push_back(*surfIterator);
+ }
+ }
+ }
+ else if (entities[ie]->geomType() == GEntity::DiscreteSurface) {
+ _ptFinalEntityList.push_back(dynamic_cast<GFace*>(entities[ie]));
+ }
+ }
#endif
- }
-
-//========================================================================================================
-
-//INITIALIZATION OF THE MAP AND RENUMBERING OF THE SELECTED ENTITIES:
+}
+// initialization of the map and renumbering of the selected entities
void Curvature::initializeMap()
{
-
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
- {
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i){
GFace* face = _ptFinalEntityList[i];
-
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
- {
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
MElement *e = face->getMeshElement(iElem);
const int E = e->getNum();
_ElementToInt[E] = 1;
-
- const int A = e->getVertex(0)->getNum(); //Pointer to 1st vertex of the triangle
+ const int A = e->getVertex(0)->getNum(); // pointer to 1st vertex of the triangle
const int B = e->getVertex(1)->getNum();
const int C = e->getVertex(2)->getNum();
-
_VertexToInt[A] = 1;
_VertexToInt[B] = 1;
_VertexToInt[C] = 1;
}
}
- /// Set up a new numbering of chosen vertices and triangles
+ // set up a new numbering of chosen vertices and triangles
int idx = 0;
// map between the pointer to vertex and the new numbering of the vertex
@@ -143,20 +110,19 @@ void Curvature::initializeMap()
// map between the pointer to element and the new numbering of the element
std::map<int,int>::iterator element_iterator;
- for (vertex_iterator = _VertexToInt.begin() ; vertex_iterator !=_VertexToInt.end() ; ++ vertex_iterator, ++idx)
+ for (vertex_iterator = _VertexToInt.begin(); vertex_iterator !=_VertexToInt.end();
+ ++ vertex_iterator, ++idx)
(*vertex_iterator).second = idx;
idx = 0;
-
- for (element_iterator = _ElementToInt.begin() ; element_iterator !=_ElementToInt.end() ; ++ element_iterator, ++idx)
+
+ for (element_iterator = _ElementToInt.begin(); element_iterator != _ElementToInt.end();
+ ++ element_iterator, ++idx)
(*element_iterator).second = idx;
}
-//========================================================================================================
-
-//COMPUTE THE NORMAL AT THE VERTEX AND THE AREA AROUND
-
+// compute the normal at the vertex and the area around
void Curvature::computeVertexNormals()
{
SVector3 vector_AB;
@@ -166,14 +132,12 @@ void Curvature::computeVertexNormals()
_VertexArea.resize(_ElementToInt.size() );
_VertexNormal.resize(_VertexToInt.size());
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
- {
+ for (int i = 0; i < _ptFinalEntityList.size(); ++i){
// face is a pointer to one surface of the group "FinalEntityList"
GFace* face = _ptFinalEntityList[i];
//Loop over the element all the element of the "myTag"-surface
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
- {
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
// Pointer to one element
MElement *e = face->getMeshElement(iElem);
// The NEW tag of the corresponding element
@@ -196,13 +160,15 @@ void Curvature::computeVertexNormals()
// Area of the triangles:
_TriangleArea[E] = 0.5*cross.norm();
- // std::cout << "The area of the triangle nr: " << e->getNum() << " is: "<< TriangleArea[E] << std::endl;
+ // std::cout << "The area of the triangle nr: " << e->getNum()
+ // << " is: "<< TriangleArea[E] << std::endl;
_VertexArea[V0] += _TriangleArea[E];
_VertexArea[V1] += _TriangleArea[E];
_VertexArea[V2] += _TriangleArea[E];
- _VertexNormal[V0] += cross; //here we are actually computing the unit normal vector per vertex
+ // here we are actually computing the unit normal vector per vertex
+ _VertexNormal[V0] += cross;
_VertexNormal[V1] += cross;
_VertexNormal[V2] += cross;
@@ -210,18 +176,14 @@ void Curvature::computeVertexNormals()
} //Loop over _ptFinalEntityList
- ///////Normalize the vertex-normals.
- for (unsigned int n = 0; n < _VertexToInt.size(); ++ n)
- {
- _VertexNormal[n].normalize();
- }
+ // Normalize the vertex-normals.
+ for (unsigned int n = 0; n < _VertexToInt.size(); ++ n){
+ _VertexNormal[n].normalize();
+ }
}
-//========================================================================================================
-
void Curvature::curvatureTensor()
{
-
STensor3 TempTensor;
STensor3 TijABTensorProduct;
SVector3 TijAB;
@@ -229,23 +191,21 @@ void Curvature::curvatureTensor()
_CurveTensor.resize(_VertexToInt.size());
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
- {
- GFace* face = _ptFinalEntityList[i]; //face is a pointer to one surface of the group "FinalEntityList"
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i){
+ // face is a pointer to one surface of the group "FinalEntityList"
+ GFace* face = _ptFinalEntityList[i];
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
- {
+ //Loop over the element all the element of the "myTag"-surface
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
MElement *e = face->getMeshElement(iElem); //Pointer to one element
const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
- for (unsigned int i = 0; i<3; ++i) // Loop over the 3 edges of each element
- {
-
- MVertex* A = e->getVertex(i); //Pointer to 1st vertex of the edge A-to-B
- MVertex* B = e->getVertex((i+1)%3); //Pointer to 2nd vertex of the edge A-to-B
-
- const int V0 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the edge A-to-B
- const int V1 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the edge A-to-B
+ // Loop over the 3 edges of each element
+ for (unsigned int i = 0; i<3; ++i){
+ MVertex* A = e->getVertex(i); //Pointer to 1st vertex of the edge A-to-B
+ MVertex* B = e->getVertex((i+1)%3); //Pointer to 2nd vertex of the edge A-to-B
+ const int V0 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the edge A-to-B
+ const int V1 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the edge A-to-B
//Weight for triangle-i-th's contribution to the shape tensor:
const double Wij0 = _TriangleArea[E] / (2 * _VertexArea[V0]);
@@ -253,7 +213,7 @@ void Curvature::curvatureTensor()
//Approximate Curvature "kappa" along some tangential vector T:
vector_AB = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z() );
- const double k_nominator0 = dot(_VertexNormal[V0], vector_AB); //Dot-product of the 2 vectors
+ const double k_nominator0 = dot(_VertexNormal[V0], vector_AB);
const double k_nominator1 = -dot(_VertexNormal[V1], vector_AB);
const double coeff = 2.0/vector_AB.normSq(); //normSq is the norm to the power 2
@@ -268,11 +228,9 @@ void Curvature::curvatureTensor()
TempTensor(1,1) += 1.0;
TempTensor(2,2) += 1.0;
- for (int m = 0; m<3; ++m)
- {
+ for (int m = 0; m<3; ++m){
TijAB(m) = 0.0;
- for (int n = 0; n<3; ++n)
- {
+ for (int n = 0; n<3; ++n){
TijAB(m) += TempTensor(m,n) * vector_AB(n);
}
}
@@ -288,11 +246,9 @@ void Curvature::curvatureTensor()
TempTensor(1,1) += 1.0;
TempTensor(2,2) += 1.0;
- for (int m = 0; m<3; ++m)
- {
+ for (int m = 0; m<3; ++m){
TijAB(m) = 0.0;
- for (int n = 0; n<3; ++n)
- {
+ for (int n = 0; n<3; ++n){
TijAB(m) += TempTensor(m,n) * vector_AB(n);
}
}
@@ -309,8 +265,6 @@ void Curvature::curvatureTensor()
}//End of method
-//========================================================================================================
-
void Curvature::computeCurvature_Simple()
{
SVector3 vector_E;
@@ -327,8 +281,8 @@ void Curvature::computeCurvature_Simple()
_VertexCurve.resize(_VertexToInt.size());
- for (unsigned int n = 0; n < _VertexToInt.size(); ++ n) //Loop over the vertex
- {
+ //Loop over the vertex
+ for (unsigned int n = 0; n < _VertexToInt.size(); ++ n){
vector_E = SVector3(1,0,0);
vector_A = vector_E + _VertexNormal[n];
vector_B = vector_E - _VertexNormal[n];
@@ -336,12 +290,10 @@ void Curvature::computeCurvature_Simple()
const double MagA = vector_A.norm();
const double MagB = vector_B.norm();
- if (MagB > MagA)
- {
+ if (MagB > MagA){
vector_Wvi = vector_B;
}
- else
- {
+ else{
vector_Wvi = vector_A;
}
@@ -354,14 +306,12 @@ void Curvature::computeCurvature_Simple()
Qvi(2,2) += 1.0;
//Transpose the matrix:
- for (int i = 0; i<3; ++i)
- {
- for (int j = 0; j<3; ++j)
- {
+ for (int i = 0; i<3; ++i){
+ for (int j = 0; j<3; ++j){
QviT(i,j) = Qvi(j,i);
- }
- }
-
+ }
+ }
+
QviT *= _CurveTensor[n];
QviT *=Qvi;
Holder = QviT;
@@ -372,8 +322,7 @@ void Curvature::computeCurvature_Simple()
const double C = Holder(1,1)*Holder(2,2) - Holder(1,2)*Holder(2,1);
const double Delta = std::sqrt(B*B-4*A*C);
- if((B*B-4.*A*C) < 0.0)
- {
+ if((B*B-4.*A*C) < 0.0){
std::cout << "WARNING: negative discriminant: " << B*B-4.*A*C << std::endl;
}
@@ -386,10 +335,6 @@ void Curvature::computeCurvature_Simple()
}
}
-//========================================================================================================
-
-//COMPUTE THE NORMAL AT THE VERTEX AND THE AREA AROUND
-
void Curvature::computeRusinkiewiczNormals()
{
SVector3 vector_AB;
@@ -399,14 +344,12 @@ void Curvature::computeRusinkiewiczNormals()
_TriangleArea.resize(_ElementToInt.size() );
_VertexNormal.resize(_VertexToInt.size());
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
- {
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i){
// face is a pointer to one surface of the group "FinalEntityList"
GFace* face = _ptFinalEntityList[i];
//Loop over the element all the element of the "myTag"-surface
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
- {
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
// Pointer to one element
MElement *e = face->getMeshElement(iElem);
const int E = _ElementToInt[e->getNum()];
@@ -416,7 +359,7 @@ void Curvature::computeRusinkiewiczNormals()
MVertex* B = e->getVertex(1);
MVertex* C = e->getVertex(2);
- const int V0 = _VertexToInt[A->getNum()]; //The new number of the vertex
+ const int V0 = _VertexToInt[A->getNum()]; //The new number of the vertex
const int V1 = _VertexToInt[B->getNum()];
const int V2 = _VertexToInt[C->getNum()];
@@ -440,18 +383,16 @@ void Curvature::computeRusinkiewiczNormals()
} // end of loop over elements of one face
} //Loop over _ptFinalEntityList
-
- ///////Normalize the vertex-normals.
- for (unsigned int n = 0; n < _VertexToInt.size(); ++ n)
- {
- _VertexNormal[n].normalize();
- }
+
+ // Normalize the vertex-normals.
+ for (unsigned int n = 0; n < _VertexToInt.size(); ++ n){
+ _VertexNormal[n].normalize();
+ }
}
-//========================================================================================================
// Compute per-vertex point areas
-void Curvature::computePointareas(){
-
+void Curvature::computePointareas()
+{
SVector3 e[3];
SVector3 l2;
SVector3 ew;
@@ -461,10 +402,11 @@ void Curvature::computePointareas(){
for (int i = 0; i< _ptFinalEntityList.size(); ++i)
{
- GFace* face = _ptFinalEntityList[i]; //face is a pointer to one surface of the group "FinalEntityList"
+ //face is a pointer to one surface of the group "FinalEntityList"
+ GFace* face = _ptFinalEntityList[i];
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
- {
+ //Loop over the element all the element of the "myTag"-surface
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
MElement *E = face->getMeshElement(iElem); //Pointer to one element
// The NEW tag of the corresponding element
const int EIdx = _ElementToInt[E->getNum()];
@@ -474,7 +416,7 @@ void Curvature::computePointareas(){
MVertex* C = E->getVertex(2);
//Edges
- e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()); //vector side of a triangilar element
+ e[0] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z());
e[1] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
e[2] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());
@@ -488,26 +430,22 @@ void Curvature::computePointareas(){
l2[1] * (l2[2] + l2[0] - l2[1]),
l2[2] * (l2[0] + l2[1] - l2[2]) );
- if (ew[0] <= 0.0)
- {
+ if (ew[0] <= 0.0){
_cornerareas[EIdx][1] = -0.25 * l2[2] * area / dot(e[0], e[2]);
_cornerareas[EIdx][2] = -0.25 * l2[1] * area / dot(e[0], e[1]);
_cornerareas[EIdx][0] = area - _cornerareas[EIdx][1] - _cornerareas[EIdx][2];
}
- else if (ew[1] <= 0.0)
- {
+ else if (ew[1] <= 0.0){
_cornerareas[EIdx][2] = -0.25 * l2[0] * area / dot(e[1], e[0]);
_cornerareas[EIdx][0] = -0.25 * l2[2] * area / dot(e[1], e[2]);
_cornerareas[EIdx][1] = area - _cornerareas[EIdx][2] - _cornerareas[EIdx][0];
}
- else if (ew[2] <= 0.0)
- {
+ else if (ew[2] <= 0.0){
_cornerareas[EIdx][0] = -0.25 * l2[1] * area / dot(e[2], e[1]);
_cornerareas[EIdx][1] = -0.25 * l2[0] * area / dot(e[2], e[0]);
_cornerareas[EIdx][2] = area - _cornerareas[EIdx][0] - _cornerareas[EIdx][1];
}
- else
- {
+ else{
float ewscale = 0.5 * area / (ew[0] + ew[1] + ew[2]);
for (int j = 0; j < 3; j++)
_cornerareas[EIdx][j] = ewscale * (ew[(j+1)%3] + ew[(j+2)%3]);
@@ -526,23 +464,23 @@ void Curvature::computePointareas(){
} //End of loop over iElem
-// std::cout << "_pointareas.size = " << _pointareas.size() << std::endl;
-// std::cout << "_cornerareas.size = " << _cornerareas.size() << std::endl;
+ // std::cout << "_pointareas.size = " << _pointareas.size() << std::endl;
+ // std::cout << "_cornerareas.size = " << _cornerareas.size() << std::endl;
} //End of loop over _ptFinalEntityList
-} //End of the method "computePointareas"
-
+}
-//========================================================================================================
//Rotate a coordinate system to be perpendicular to the given normal
-void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v, const SVector3 &new_norm, SVector3 &new_u, SVector3 &new_v){
-
+void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v,
+ const SVector3 &new_norm, SVector3 &new_u,
+ SVector3 &new_v)
+{
new_u = old_u;
new_v = old_v;
SVector3 old_norm = crossprod(old_u, old_v);
double ndot = dot(old_norm, new_norm);
-// if (unlikely(ndot <= -1.0f))
+// if (unlikely(ndot <= -1.0f))
if (ndot <= -1.0f)
{
new_u = -1.0*new_u;
@@ -556,16 +494,14 @@ void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v, cons
new_v -= dperp*dot(new_v, perp_old);
}
-//========================================================================================================
-
-//Project a curvature tensor from the basis spanned by old_u and old_v
-//(which are assumed to be unit-length and perpendicular) to the new_u
-//and new_v basis
-
-void Curvature::proj_curv( const SVector3 &old_u, const SVector3 &old_v,
+// Project a curvature tensor from the basis spanned by old_u and
+// old_v (which are assumed to be unit-length and perpendicular) to
+// the new_u and new_v basis
+void Curvature::proj_curv(const SVector3 &old_u, const SVector3 &old_v,
double old_ku, double old_kuv, double old_kv,
const SVector3 &new_u, const SVector3 &new_v,
- double &new_ku, double &new_kuv, double &new_kv){
+ double &new_ku, double &new_kuv, double &new_kv)
+{
SVector3 r_new_u;
SVector3 r_new_v;
rot_coord_sys(new_u, new_v, crossprod(old_u,old_v), r_new_u, r_new_v);
@@ -580,16 +516,14 @@ void Curvature::proj_curv( const SVector3 &old_u, const SVector3 &old_v,
new_kv = old_ku*u2*u2 + old_kuv*(2.0f * u2*v2) + old_kv*v2*v2;
}
-
-//========================================================================================================
-
-//Given a curvature tensor, find principal directions and curvatures
-//Makes sure that pdir1 and pdir2 are perpendicular to normal
-
+// Given a curvature tensor, find principal directions and curvatures
+// Makes sure that pdir1 and pdir2 are perpendicular to normal
void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
- double ku, double kuv, double kv,
- const SVector3 &new_norm,
- SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2){
+ double ku, double kuv, double kv,
+ const SVector3 &new_norm,
+ SVector3 &pdir1, SVector3 &pdir2,
+ double &k1, double &k2)
+{
SVector3 r_old_u;
SVector3 r_old_v;
@@ -599,7 +533,7 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
rot_coord_sys(old_u, old_v, new_norm, r_old_u, r_old_v);
-// if(unlikely(kuv !=0.0f))
+// if(unlikely(kuv !=0.0f))
if(kuv !=0.0) {
//Jacobi rotation to diagonalize
double h= 0.5*(kv -ku)/ kuv;
@@ -613,20 +547,18 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
k1 = ku - tt *kuv;
k2 = kv + tt *kuv;
- if(std::abs(k1) >= std::abs(k2)) {
+ if(std::abs(k1) >= std::abs(k2)) {
pdir1 = c*r_old_u - s*r_old_v;
}
- else {
+ else {
std::swap(k1,k2);
pdir1 = s*r_old_u + c*r_old_v;
}
pdir2 = crossprod(new_norm, pdir1);
}
-//========================================================================================================
void Curvature::computeCurvature(GModel* model, typeOfCurvature typ)
{
-
_model = model;
double t0 = Cpu();
@@ -643,7 +575,6 @@ void Curvature::computeCurvature(GModel* model, typeOfCurvature typ)
writeToPosFile("curvature.pos");
writeToVtkFile("curvature.vtk");
-
}
void Curvature::computeCurvature_Rusinkiewicz(int isMax)
@@ -681,44 +612,41 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
for (int i = 0; i< _ptFinalEntityList.size(); ++i)
{
- GFace* face = _ptFinalEntityList[i]; //face is a pointer to one surface of the group "FinalEntityList"
+ //face is a pointer to one surface of the group "FinalEntityList"
+ GFace* face = _ptFinalEntityList[i];
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
- {
- MElement *E = face->getMeshElement(iElem); //Pointer to one element
+ //Loop over the element all the element of the "myTag"-surface
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ MElement *E = face->getMeshElement(iElem); // Pointer to one element
- MVertex* A = E->getVertex(0); //Pointers to vertices of triangle
+ MVertex* A = E->getVertex(0); // Pointers to vertices of triangle
MVertex* B = E->getVertex(1);
MVertex* C = E->getVertex(2);
- const int V0 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
+ const int V0 = _VertexToInt[A->getNum()]; // Tag of the 1st vertex of the triangle
const int V1 = _VertexToInt[B->getNum()];
const int V2 = _VertexToInt[C->getNum()];
- ///Set up an initial coordinate system per vertex:
-
+ //Set up an initial coordinate system per vertex:
_pdir1[V0] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());
_pdir1[V1] = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z());
_pdir1[V2] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
}
}
- for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
- {
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex){
_pdir1[ivertex] = crossprod(_pdir1[ivertex], _VertexNormal[ivertex]);
_pdir1[ivertex].normalize();
_pdir2[ivertex] = crossprod(_VertexNormal[ivertex], _pdir1[ivertex]);
}
// Compute curvature per face:
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
- {
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i){
//face is a pointer to one surface of the group "FinalEntityList"
GFace* face = _ptFinalEntityList[i];
//Loop over all elements of this face:
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
- {
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
MElement *E = face->getMeshElement(iElem); //Pointer to one element
// The NEW tag of the corresponding element
const int EIdx = _ElementToInt[E->getNum()];
@@ -731,31 +659,31 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
e[1] = SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z());
e[2] = SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z());
- //SVector3 e[3] = {SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z()), SVector3(A->x() - C->x(), A->y() - C->y(), A->z() - C->z()), SVector3(B->x() - A->x(), B->y() - A->y(), B->z() - A->z()) };
+ // SVector3 e[3] = {SVector3(C->x() - B->x(), C->y() - B->y(),
+ // C->z() - B->z()), SVector3(A->x() - C->x(), A->y() - C->y(),
+ // A->z() - C->z()), SVector3(B->x() - A->x(), B->y() - A->y(),
+ // B->z() - A->z()) };
- //N-T-B coordinate system per face
+ // N-T-B coordinate system per face
t = e[0];
t.normalize();
n = crossprod( e[0], e[1]);
b = crossprod(n, t);
b.normalize();
- //Estimate curvature based on variations of normals along edges:
- //intialization:
+ // Estimate curvature based on variations of normals along
+ // edges: intialization:
m = SVector3(0.0, 0.0, 0.0);
//maybe double m[3] = { 0.0, 0.0, 0.0 };
// w *= 0.0; //Reset w to zero
- for (int i = 0; i< 3; ++i)
- {
- for (int j = 0; j<3; ++j)
- {
+ for (int i = 0; i< 3; ++i){
+ for (int j = 0; j<3; ++j){
w(i,j) = 0.0;
}
}
//filling:
- for (int j = 0; j< 3; ++j)
- {
+ for (int j = 0; j< 3; ++j){
u = dot(e[j], t);
v = dot(e[j], b);
@@ -782,23 +710,21 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
w(1,1) = w(0,0) + w(2,2);
w(1,2) = w(0,1);
- //Least Squares Solution
+ // Least Squares Solution
double diag[3];
- if (!ldltdc(w, diag))
- {
+ if (!ldltdc(w, diag)){
std::cout << "ldltdc failed" << std::endl;
continue;
}
ldltsl(w, diag, m, m);
- //Push it back out to the vertices
- for (int j = 0; j< 3; ++j)
- {
+ // Push it back out to the vertices
+ for (int j = 0; j< 3; ++j){
const int old_vj = E->getVertex(j)->getNum();
const int vj = _VertexToInt[old_vj];
proj_curv(t, b, m[0], m[1], m[2], _pdir1[vj], _pdir2[vj], c1, c12, c2);
wt = _cornerareas[EIdx][j]/_pointareas[vj];
-// wt = 1.0;
+ // wt = 1.0;
_curv1[vj] += wt*c1;
_curv12[vj] += wt*c12;
@@ -809,17 +735,17 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
} //End of loop over "_ptFinalEntityList"
-
//Compute principal directions and curvatures at each vertex
- for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex) {
- diagonalize_curv(_pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex], _curv12[ivertex], _curv2[ivertex],
- _VertexNormal[ivertex], _pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex], _curv2[ivertex]);
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex) {
+ diagonalize_curv(_pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex],
+ _curv12[ivertex], _curv2[ivertex],
+ _VertexNormal[ivertex], _pdir1[ivertex], _pdir2[ivertex],
+ _curv1[ivertex], _curv2[ivertex]);
}
-
- _VertexCurve.resize( _VertexToInt.size() );
+
+ _VertexCurve.resize(_VertexToInt.size());
for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex){
-
if (isMax){
_VertexCurve[ivertex] = std::max(fabs(_curv1[ivertex]), fabs(_curv2[ivertex]));
}
@@ -833,15 +759,14 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax)
}
_alreadyComputedCurvature = true;
-} //End of the "computeCurvature_Rusinkiewicz" method
-
+}
void Curvature::computeCurvature_RBF()
{
retrieveCompounds();
initializeMap();
- //fill set of MVertex
+ // fill set of MVertex
std::set<MVertex*> allNodes;
for (int i = 0; i< _ptFinalEntityList.size(); ++i) {
GFaceCompound* face = (GFaceCompound*)_ptFinalEntityList[i];
@@ -853,7 +778,7 @@ void Curvature::computeCurvature_RBF()
}
}
- //bounding box
+ // bounding box
SBoundingBox3d bb;
std::vector<SPoint3> vertices;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
@@ -863,18 +788,18 @@ void Curvature::computeCurvature_RBF()
}
double sizeBox = norm(SVector3(bb.max(), bb.min()));
- //compure curvature RBF
+ // compure curvature RBF
std::map<MVertex*, SVector3> _normals;
std::vector<MVertex*> _ordered;
std::map<MVertex*, double> curvRBF;
- //GLOBAL
+ // GLOBAL
GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered);
_rbf->computeCurvature(_rbf->getXYZ(),curvRBF);
- //LOCAL FD
- //GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered, true);
- //_rbf->computeLocalCurvature(_rbf->getXYZ(),curvRBF);
+ // LOCAL FD
+ // GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered, true);
+ // _rbf->computeLocalCurvature(_rbf->getXYZ(),curvRBF);
- //fill vertex curve
+ // fill vertex curve
_VertexCurve.resize( _VertexToInt.size() );
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
@@ -883,57 +808,58 @@ void Curvature::computeCurvature_RBF()
vertexIterator = _VertexToInt.find( v->getNum() );
if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
_VertexCurve[V0] = curvRBF[v];
- }
+ }
_alreadyComputedCurvature = true;
-} //End of the "computeCurvature_RBF" method
-
-
- //========================================================================================================
-
-void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs)
- {
- MVertex* A = triangle->getVertex(0);
- MVertex* B = triangle->getVertex(1);
- MVertex* C = triangle->getVertex(2);
-
- int V0 = 0;
- int V1 = 0;
- int V2 = 0;
-
- std::map<int,int>::iterator vertexIterator;
- vertexIterator = _VertexToInt.find( A->getNum() );
- if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
- else
- std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
-
- vertexIterator = _VertexToInt.find( B->getNum() );
- if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
- else
- std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
-
- vertexIterator = _VertexToInt.find( C->getNum() );
- if ( vertexIterator != _VertexToInt.end() ) V2 = (*vertexIterator).second;
- else
- std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl;
-
- if (isAbs){
- c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
- c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
- c2 = std::abs(_VertexCurve[V2]); //Mean curvature in vertex 2
- }
- else{
- c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
- c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
- c2 = _VertexCurve[V2]; //Mean curvature in vertex 2
- }
+}
+void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1,
+ double& c2, int isAbs)
+{
+ MVertex* A = triangle->getVertex(0);
+ MVertex* B = triangle->getVertex(1);
+ MVertex* C = triangle->getVertex(2);
+
+ int V0 = 0;
+ int V1 = 0;
+ int V2 = 0;
+
+ std::map<int,int>::iterator vertexIterator;
+ vertexIterator = _VertexToInt.find( A->getNum() );
+ if (vertexIterator != _VertexToInt.end()) V0 = (*vertexIterator).second;
+ else
+ std::cout << "Didn't find vertex with number " << A->getNum()
+ << " in _VertextToInt !" << std::endl;
+
+ vertexIterator = _VertexToInt.find( B->getNum() );
+ if (vertexIterator != _VertexToInt.end()) V1 = (*vertexIterator).second;
+ else
+ std::cout << "Didn't find vertex with number " << B->getNum()
+ << " in _VertextToInt !" << std::endl;
+
+ vertexIterator = _VertexToInt.find( C->getNum() );
+ if ( vertexIterator != _VertexToInt.end() ) V2 = (*vertexIterator).second;
+ else
+ std::cout << "Didn't find vertex with number " << C->getNum()
+ << " in _VertextToInt !" << std::endl;
+
+ if (isAbs){
+ c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
+ c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
+ c2 = std::abs(_VertexCurve[V2]); //Mean curvature in vertex 2
}
+ else{
+ c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
+ c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
+ c2 = _VertexCurve[V2]; //Mean curvature in vertex 2
+ }
+
+}
-//========================================================================================================
-
-void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs)
+void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax,
+ SVector3* dMin, double* cMax,
+ double* cMin, int isAbs)
{
MVertex* A = triangle->getVertex(0);
MVertex* B = triangle->getVertex(1);
@@ -945,19 +871,22 @@ void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3*
std::map<int,int>::iterator vertexIterator;
vertexIterator = _VertexToInt.find( A->getNum() );
- if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
+ if (vertexIterator != _VertexToInt.end()) V0 = (*vertexIterator).second;
else
- std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
+ std::cout << "Didn't find vertex with number " << A->getNum()
+ << " in _VertextToInt !" << std::endl;
vertexIterator = _VertexToInt.find( B->getNum() );
- if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
+ if (vertexIterator != _VertexToInt.end()) V1 = (*vertexIterator).second;
else
- std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
+ std::cout << "Didn't find vertex with number " << B->getNum()
+ << " in _VertextToInt !" << std::endl;
vertexIterator = _VertexToInt.find( C->getNum() );
- if ( vertexIterator != _VertexToInt.end() ) V2 = (*vertexIterator).second;
+ if (vertexIterator != _VertexToInt.end()) V2 = (*vertexIterator).second;
else
- std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl;
+ std::cout << "Didn't find vertex with number " << C->getNum()
+ << " in _VertextToInt !" << std::endl;
if (isAbs){
dMax[0] = _pdir1[V0];
@@ -975,10 +904,8 @@ void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3*
cMin[0] = std::abs(_curv2[V0]);
cMin[1] = std::abs(_curv2[V1]);
cMin[2] = std::abs(_curv2[V2]);
-
}
else{
-
dMax[0] = _pdir1[V0];
dMax[1] = _pdir1[V1];
dMax[2] = _pdir1[V2];
@@ -997,45 +924,42 @@ void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3*
}
}
-
-
- //========================================================================================================
-
void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs)
- {
- MVertex* A = edge->getVertex(0);
- MVertex* B = edge->getVertex(1);
-
- int V0 = 0;
- int V1 = 0;
-
- std::map<int,int>::iterator vertexIterator;
-
- vertexIterator = _VertexToInt.find( A->getNum() );
- if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
- else std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
-
- vertexIterator = _VertexToInt.find( B->getNum() );
- if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
- else std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
+{
+ MVertex* A = edge->getVertex(0);
+ MVertex* B = edge->getVertex(1);
+
+ int V0 = 0;
+ int V1 = 0;
+
+ std::map<int,int>::iterator vertexIterator;
+
+ vertexIterator = _VertexToInt.find( A->getNum() );
+ if (vertexIterator != _VertexToInt.end()) V0 = (*vertexIterator).second;
+ else std::cout << "Didn't find vertex with number " << A->getNum()
+ << " in _VertextToInt !" << std::endl;
- if (isAbs){
- c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
- c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
- }
- else{
- c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
- c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
- }
-
- }
-
-//========================================================================================================
+ vertexIterator = _VertexToInt.find( B->getNum() );
+ if (vertexIterator != _VertexToInt.end()) V1 = (*vertexIterator).second;
+ else std::cout << "Didn't find vertex with number " << B->getNum()
+ << " in _VertextToInt !" << std::endl;
+
+ if (isAbs){
+ c0 = std::abs(_VertexCurve[V0]); //Mean curvature in vertex 0
+ c1 = std::abs(_VertexCurve[V1]); //Mean curvature in vertex 1
+ }
+ else{
+ c0 = _VertexCurve[V0]; //Mean curvature in vertex 0
+ c1 = _VertexCurve[V1]; //Mean curvature in vertex 1
+ }
+}
-double Curvature::getAtVertex(const MVertex *v) const {
+double Curvature::getAtVertex(const MVertex *v) const
+{
std::map<int,int>::const_iterator it = _VertexToInt.find(v->getNum());
if (it == _VertexToInt.end()) {
- Msg::Error("curvature has not been computed for vertex %i (%i)", v->getNum(), _VertexToInt.size());
+ Msg::Error("curvature has not been computed for vertex %i (%i)",
+ v->getNum(), _VertexToInt.size());
return 1;
}
return _VertexCurve[it->second];
@@ -1047,7 +971,7 @@ void Curvature::writeToPosFile( const std::string & filename)
outfile.precision(18);
outfile.open(filename.c_str());
outfile << "View \"Curvature \"{" << std::endl;
-
+
int idxelem = 0;
for (int i = 0; i< _ptFinalEntityList.size(); ++i) {
@@ -1062,12 +986,11 @@ void Curvature::writeToPosFile( const std::string & filename)
MVertex* B = e->getVertex(1);
MVertex* C = e->getVertex(2);
- const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
- const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
- const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
-
- //Here is printing the triplet X-Y-Z of each vertex:
- //*************************************************
+ const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
+ const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
+ const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
+
+ // Here is printing the triplet X-Y-Z of each vertex:
outfile << "ST("; //VT = vector triangles //ST = scalar triangle
outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
@@ -1077,33 +1000,31 @@ void Curvature::writeToPosFile( const std::string & filename)
outfile << ")";
outfile <<"{";
- //Here is printing the 3 components of the normal vector for each vertex:
- //**********************************************************************
-
-// outfile << VertexNormal[V1].x() << ","<< VertexNormal[V1].y() << ","<< VertexNormal[V1].z() << ",";
-// outfile << VertexNormal[V2].x() << ","<< VertexNormal[V2].y() << ","<< VertexNormal[V2].z() << ",";
-// outfile << VertexNormal[V3].x() << ","<< VertexNormal[V3].y() << ","<< VertexNormal[V3].z();
+ // Here is printing the 3 components of the normal vector for each vertex:
+ //outfile << VertexNormal[V1].x() << ","<< VertexNormal[V1].y() << ","
+ // << VertexNormal[V1].z() << ",";
+ //outfile << VertexNormal[V2].x() << ","<< VertexNormal[V2].y() << ","
+ // << VertexNormal[V2].z() << ",";
+ //outfile << VertexNormal[V3].x() << ","<< VertexNormal[V3].y() << ","
+ // << VertexNormal[V3].z();
+
outfile << _VertexCurve[V1] << "," << _VertexCurve[V2] << "," << _VertexCurve[V3];
outfile << "};" << std::endl;
idxelem++;
- } //Loop over elements
-
-} // Loop over ptFinalEntityList
-
-outfile << "};" << std::endl;
+ }
+ }
-outfile.close();
+ outfile << "};" << std::endl;
+
+ outfile.close();
}
-//========================================================================================================
-
void Curvature::writeToVtkFile( const std::string & filename)
{
-
std::ofstream outfile;
outfile.precision(18);
outfile.open(filename.c_str());
@@ -1116,22 +1037,20 @@ void Curvature::writeToVtkFile( const std::string & filename)
outfile << "POINTS " << npoints << " double" << std::endl;
- /// Build a table of coordinates
- /// Loop over all elements and look at the 'old' (not necessarily continuous) numbers of vertices
- /// Get the 'new' index of each vertex through _VertexToInt and the [x,y,z] coordinates of this vertex
- /// Store them in coordx,coordy and coordz
-
+ // Build a table of coordinates
+ // - Loop over all elements and look at the 'old' (not necessarily
+ // continuous) numbers of vertices
+ // - Get the 'new' index of each vertex through _VertexToInt and the
+ // [x,y,z] coordinates of this vertex
+ // - Store them in coordx,coordy and coordz
std::vector<VtkPoint> coord;
-
coord.resize(npoints);
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
- {
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i){
GFace* face = _ptFinalEntityList[i];
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
- {
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
MElement *e = face->getMeshElement(iElem);
MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
@@ -1160,25 +1079,23 @@ void Curvature::writeToVtkFile( const std::string & filename)
}
}
- for (int v = 0; v < npoints; ++v)
- {
+ for (int v = 0; v < npoints; ++v){
outfile << coord[v].x << " " << coord[v].y << " " << coord[v].z << std::endl;
}
- /// Empty the array 'coord' to free the memory
- /// Point coordinates will not be needed anymore
+ // Empty the array 'coord' to free the memory
+ // Point coordinates will not be needed anymore
coord.clear();
- /// Write the cell connectivity
+ // Write the cell connectivity
- outfile << std::endl << "CELLS " << _ElementToInt.size() << " " << 4*_ElementToInt.size() << std::endl;
+ outfile << std::endl << "CELLS " << _ElementToInt.size() << " "
+ << 4*_ElementToInt.size() << std::endl;
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
- {
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i){
GFace* face = _ptFinalEntityList[i];
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
- {
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
MElement *e = face->getMeshElement(iElem);
MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
@@ -1196,32 +1113,25 @@ void Curvature::writeToVtkFile( const std::string & filename)
outfile << "3 " << newIdxA << " " << newIdxB << " " << newIdxC << std::endl;
}
}
-
+
outfile << std::endl << "CELL_TYPES " << _ElementToInt.size() << std::endl;
- for(int ie = 0; ie < _ElementToInt.size(); ++ie)
- {
+ for(int ie = 0; ie < _ElementToInt.size(); ++ie){
outfile << "5" << std::endl; //Triangle is element type 5 in vtk
-
}
- /// Write the curvature values as vtk 'point data'
+ // Write the curvature values as vtk 'point data'
outfile << std::endl << "POINT_DATA " << npoints << std::endl;
outfile << "SCALARS curvature float 1" << std::endl;
outfile << "LOOKUP_TABLE default" << std::endl;
- for (int iv = 0; iv < npoints; ++iv)
- {
+ for (int iv = 0; iv < npoints; ++iv){
outfile << _VertexCurve[iv] << std::endl;
}
-
+
outfile.close();
-
-
}
-//========================================================================================================
-
void Curvature::writeDirectionsToPosFile( const std::string & filename)
{
std::ofstream outfile;
@@ -1229,28 +1139,26 @@ void Curvature::writeDirectionsToPosFile( const std::string & filename)
outfile.open(filename.c_str());
outfile << "View \"Curvature_DirMax \"{" << std::endl;
- for (int i = 0; i< _ptFinalEntityList.size(); ++i)
- {
- GFace* face = _ptFinalEntityList[i]; //face is a pointer to one surface of the group "FinalEntityList"
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i){
+ //face is a pointer to one surface of the group "FinalEntityList"
+ GFace* face = _ptFinalEntityList[i];
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
- {
+ //Loop over the element all the element of the "myTag"-surface
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
MElement *e = face->getMeshElement(iElem); //Pointer to one element
const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
-
+
//std::cout << "We are now looking at element Nr: " << E << std::endl;
-
+
MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
MVertex* B = e->getVertex(1);
MVertex* C = e->getVertex(2);
- const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
- const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
- const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
+ const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
+ const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
+ const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
//Here is printing the triplet X-Y-Z of each vertex:
- //*************************************************
-
outfile << "VT("; //VT = vector triangles //ST = scalar triangle
outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
@@ -1260,75 +1168,54 @@ void Curvature::writeDirectionsToPosFile( const std::string & filename)
outfile <<"{";
//Here is printing the 3 components of the curvature max direction for each vertex:
- //*********************************************************************************
-
- outfile << _pdir1[V1].x() << ","<< _pdir1[V1].y() << ","<< _pdir1[V1].z() << ",";
- outfile << _pdir1[V2].x() << ","<< _pdir1[V2].y() << ","<< _pdir1[V2].z() << ",";
- outfile << _pdir1[V3].x() << ","<< _pdir1[V3].y() << ","<< _pdir1[V3].z();
-
+ outfile << _pdir1[V1].x() << ","<< _pdir1[V1].y() << ","<< _pdir1[V1].z() << ",";
+ outfile << _pdir1[V2].x() << ","<< _pdir1[V2].y() << ","<< _pdir1[V2].z() << ",";
+ outfile << _pdir1[V3].x() << ","<< _pdir1[V3].y() << ","<< _pdir1[V3].z();
outfile << "};" << std::endl;
-
- } //Loop over elements
-
-} // Loop over ptFinalEntityList
-
-outfile << "};" << std::endl;
-
-
-//----------------------------------------------------------------------------------------------
-
-outfile << "View \"Curvature_DirMin \"{" << std::endl;
-
-for (int i = 0; i< _ptFinalEntityList.size(); ++i)
-{
- GFace* face = _ptFinalEntityList[i]; //face is a pointer to one surface of the group "FinalEntityList"
-
- for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
- {
- MElement *e = face->getMeshElement(iElem); //Pointer to one element
- const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
-
- //std::cout << "We are now looking at element Nr: " << E << std::endl;
-
- MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
- MVertex* B = e->getVertex(1);
- MVertex* C = e->getVertex(2);
-
- const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
- const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
- const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
-
- //Here is printing the triplet X-Y-Z of each vertex:
- //*************************************************
-
- outfile << "VT("; //VT = vector triangles //ST = scalar triangle
- outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
- outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
- outfile << C->x() << ","<< C->y() << "," << C->z();
-
- outfile << ")";
- outfile <<"{";
-
- //Here is printing the 3 components of the curvature min direction for each vertex:
- //*********************************************************************************
-
- outfile << _pdir2[V1].x() << ","<< _pdir2[V1].y() << ","<< _pdir2[V1].z() << ",";
- outfile << _pdir2[V2].x() << ","<< _pdir2[V2].y() << ","<< _pdir2[V2].z() << ",";
- outfile << _pdir2[V3].x() << ","<< _pdir2[V3].y() << ","<< _pdir2[V3].z();
-
-
- outfile << "};" << std::endl;
-
-} //Loop over elements
-
-} // Loop over ptFinalEntityList
-
-outfile << "};" << std::endl;
-
-outfile.close();
+ }
+ }
+
+ outfile << "};" << std::endl;
+
+ outfile << "View \"Curvature_DirMin \"{" << std::endl;
+
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i){
+ GFace* face = _ptFinalEntityList[i];
+
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ MElement *e = face->getMeshElement(iElem); //Pointer to one element
+ const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
+
+ //std::cout << "We are now looking at element Nr: " << E << std::endl;
+
+ MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
+ MVertex* B = e->getVertex(1);
+ MVertex* C = e->getVertex(2);
+
+ const int V1 = _VertexToInt[A->getNum()];
+ const int V2 = _VertexToInt[B->getNum()];
+ const int V3 = _VertexToInt[C->getNum()];
+
+ //Here is printing the triplet X-Y-Z of each vertex:
+ outfile << "VT("; //VT = vector triangles //ST = scalar triangle
+ outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
+ outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
+ outfile << C->x() << ","<< C->y() << "," << C->z();
+
+ outfile << ")";
+ outfile <<"{";
+
+ //Here is printing the 3 components of the curvature min direction for each vertex:
+ outfile << _pdir2[V1].x() << ","<< _pdir2[V1].y() << ","<< _pdir2[V1].z() << ",";
+ outfile << _pdir2[V2].x() << ","<< _pdir2[V2].y() << ","<< _pdir2[V2].z() << ",";
+ outfile << _pdir2[V3].x() << ","<< _pdir2[V3].y() << ","<< _pdir2[V3].z();
+
+ outfile << "};" << std::endl;
+ }
+ }
+
+ outfile << "};" << std::endl;
+
+ outfile.close();
}
-//========================================================================================================
-
-
-
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index 206f370cef15db9a9d7fd05b38c1e3c220c955f8..c34a909ce87d476f31e5121f02cc0c63a07ec547 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -3,190 +3,160 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
-#ifndef _CURVATUREL_H_
+#ifndef _CURVATURE_H_
#define _CURVATURE_H_
+#include <map>
+#include <vector>
#include "GModel.h"
#include "STensor3.h"
-#include<map>
-#include<vector>
-
class Curvature {
+ private:
+ typedef std::vector<GFace*> GFaceList;
+ typedef std::map<int, GFaceList >::iterator GFaceIter;
+
+ // helper type for writing VTK files
+ struct VtkPoint
+ {
+ double x;
+ double y;
+ double z;
+ };
+
+ // static member variable to implement the class curvature as a Singleton
+ static Curvature *_instance;
+ static bool _destroyed;
-private:
- typedef std::vector<GFace*> GFaceList;
- //typedef std::map<int, GEntityList >::iterator GEntityIter;
- typedef std::map<int, GFaceList >::iterator GFaceIter;
-
- //-----------------------------------------
- // HELPER TYPE FOR WRITING TO VTK FILES
- struct VtkPoint
- {
- double x;
- double y;
- double z;
- };
-
- //-----------------------------------------
- // MEMBER VARIABLES
-
- //Static member variable to implement the class curvature as a Singleton
- static Curvature *_instance;
- static bool _destroyed;
-
- //Boolean to check if the curvature has already been computed
- static bool _alreadyComputedCurvature;
-
- //Map definition
- std::map<int, int> _VertexToInt;
- std::map<int, int> _ElementToInt;
-
- //Model and list of selected entities with give physical tag:
- GModel* _model;
- GFaceList _ptFinalEntityList;
-
- //Averaged vertex normals
- std::vector<SVector3> _VertexNormal;
-
- // Vector of principal dircections
- std::vector<SVector3> _pdir1;
- std::vector<SVector3> _pdir2;
-
- // Vector of principal curvature
- std::vector<double> _curv1;
- std::vector<double> _curv2;
- std::vector<double> _curv12;
-
- // Area around point
- std::vector<double> _pointareas;
- std::vector<SVector3> _cornerareas;
-
- //Curvature Tensor per mesh vertex
- std::vector<STensor3> _CurveTensor;
-
- //Area of a triangle element:
- std::vector<double> _TriangleArea;
-
- //Area around a mesh vertex:
- std::vector<double> _VertexArea;
- std::vector<double> _VertexCurve;
+ // boolean to check if the curvature has already been computed
+ static bool _alreadyComputedCurvature;
- //-----------------------------------------
- // PRIVATE METHODS
+ // map definition
+ std::map<int, int> _VertexToInt;
+ std::map<int, int> _ElementToInt;
- //Constructor and destructor
- Curvature();
- ~Curvature();
+ // model and list of selected entities with give physical tag:
+ GModel* _model;
+ GFaceList _ptFinalEntityList;
- static void create();
- static void onDeadReference();
+ // averaged vertex normals
+ std::vector<SVector3> _VertexNormal;
- void initializeMap();
- void computeVertexNormals();
- void curvatureTensor();
- static void rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v,
- const SVector3 &new_norm, SVector3 &new_u, SVector3 &new_v);
- void proj_curv( const SVector3 &old_u, const SVector3 &old_v, double old_ku, double old_kuv,
- double old_kv, const SVector3 &new_u, const SVector3 &new_v,
- double &new_ku, double &new_kuv, double &new_kv);
- void diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
- double ku, double kuv, double kv,
- const SVector3 &new_norm,
- SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2);
- void computePointareas();
- void computeRusinkiewiczNormals();
+ // vector of principal dircections
+ std::vector<SVector3> _pdir1;
+ std::vector<SVector3> _pdir2;
- // Perform LDL^T decomposition of a symmetric positive definite matrix.
- // Like Cholesky, but no square roots. Overwrites lower triangle of matrix.
- static inline bool ldltdc(STensor3& A, double rdiag[3])
- {
- double v[2];
- for (int i = 0; i < 3; i++)
- {
+ // vector of principal curvature
+ std::vector<double> _curv1;
+ std::vector<double> _curv2;
+ std::vector<double> _curv12;
+
+ // area around point
+ std::vector<double> _pointareas;
+ std::vector<SVector3> _cornerareas;
+
+ // curvature Tensor per mesh vertex
+ std::vector<STensor3> _CurveTensor;
+
+ // area of a triangle element:
+ std::vector<double> _TriangleArea;
+
+ // area around a mesh vertex:
+ std::vector<double> _VertexArea;
+ std::vector<double> _VertexCurve;
+
+ // constructor and destructor
+ Curvature(){}
+ ~Curvature();
+ static void create();
+ static void onDeadReference();
+ void initializeMap();
+ void computeVertexNormals();
+ void curvatureTensor();
+ static void rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v,
+ const SVector3 &new_norm, SVector3 &new_u,
+ SVector3 &new_v);
+ void proj_curv(const SVector3 &old_u, const SVector3 &old_v, double old_ku,
+ double old_kuv, double old_kv, const SVector3 &new_u,
+ const SVector3 &new_v, double &new_ku, double &new_kuv,
+ double &new_kv);
+ void diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
+ double ku, double kuv, double kv,
+ const SVector3 &new_norm,
+ SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2);
+ void computePointareas();
+ void computeRusinkiewiczNormals();
+
+ // Perform LDL^T decomposition of a symmetric positive definite matrix.
+ // Like Cholesky, but no square roots. Overwrites lower triangle of matrix.
+ static inline bool ldltdc(STensor3& A, double rdiag[3])
+ {
+ double v[2];
+ for (int i = 0; i < 3; i++){
+ for (int k = 0; k < i; k++)
+ v[k] = A(i,k) * rdiag[k];
+ for (int j = i; j < 3; j++){
+ double sum = A(i,j);
for (int k = 0; k < i; k++)
- v[k] = A(i,k) * rdiag[k];
- for (int j = i; j < 3; j++)
- {
- double sum = A(i,j);
- for (int k = 0; k < i; k++)
- sum -= v[k] * A(j,k);
- if (i == j)
- {
- //if (unlikely(sum <= T(0)))
- if (sum <= 0.0)
- return false;
- rdiag[i] = 1.0 / sum;
- }
- else
- {
- A(j,i) = sum;
- }
+ sum -= v[k] * A(j,k);
+ if (i == j){
+ //if (unlikely(sum <= T(0)))
+ if (sum <= 0.0)
+ return false;
+ rdiag[i] = 1.0 / sum;
+ }
+ else{
+ A(j,i) = sum;
}
}
-
- return true;
}
- // Solve Ax=B after ldltdc
- static inline void ldltsl(STensor3& A, double rdiag[3], double B[3], double x[3])
- {
- int i;
- for (i = 0; i < 3; i++)
- {
- double sum = B[i];
- for (int k = 0; k < i; k++)
- sum -= A(i,k) * x[k];
- x[i] = sum * rdiag[i];
- }
- for (i = 3 - 1; i >= 0; i--)
- {
- double sum = 0;
- for (int k = i + 1; k < 3; k++)
- sum += A(k,i) * x[k];
- x[i] -= sum * rdiag[i];
- }
+ return true;
+ }
+ // Solve Ax=B after ldltdc
+ static inline void ldltsl(STensor3& A, double rdiag[3], double B[3], double x[3])
+ {
+ int i;
+ for (i = 0; i < 3; i++){
+ double sum = B[i];
+ for (int k = 0; k < i; k++)
+ sum -= A(i,k) * x[k];
+ x[i] = sum * rdiag[i];
}
-
-public:
-
- typedef enum {RUSIN=1,RBF=2, SIMPLE=3} typeOfCurvature;
+ for (i = 3 - 1; i >= 0; i--){
+ double sum = 0;
+ for (int k = i + 1; k < 3; k++)
+ sum += A(k,i) * x[k];
+ x[i] -= sum * rdiag[i];
+ }
+ }
+
+ public:
+ typedef enum {RUSIN=1, RBF=2, SIMPLE=3} typeOfCurvature;
static Curvature& getInstance();
static bool valueAlreadyComputed();
-
void retrieveCompounds();
double getAtVertex(const MVertex *v) const;
- //void retrievePhysicalSurfaces(const std::string & face_tag);
-
void computeCurvature(GModel* model, typeOfCurvature typ);
-
/// The following function implements algorithm from:
/// Implementation of an Algorithm for Approximating the Curvature Tensor
/// on a Triangular Surface Mesh in the Vish Environment
/// Edwin Matthews, Werner Benger, Marcel Ritter
void computeCurvature_Simple();
-
/// The following function implements algorithm from:
/// Estimating Curvatures and Their Derivatives on Triangle Meshes
/// Szymon Rusinkiewicz, Princeton University
/// Code taken from Rusinkiewicz' 'trimesh2' library
void computeCurvature_Rusinkiewicz(int isMax=0);
-
void computeCurvature_RBF();
-
- void triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs=0);
- void triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs=0);
-
+ void triangleNodalValues(MTriangle* triangle, double& c0, double& c1,
+ double& c2, int isAbs=0);
+ void triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax,
+ SVector3* dMin, double* cMax, double* cMin,
+ int isAbs=0);
void edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs=0);
-
void writeToPosFile( const std::string & filename);
-
void writeToVtkFile( const std::string & filename);
-
void writeDirectionsToPosFile( const std::string & filename);
-
-
-
};
-
#endif
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 6d99b8f8d1081c04a1511f67dfd725d765b475ed..9b67d26bcea6907125c7bd46bfb5cca3f3c051ea 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -78,7 +78,7 @@ void GEdgeCompound::orderEdges()
else tempv.erase(it2);
}
- //find the first GEdge and erase it from the list edges
+ //find the first GEdge and erase it from the list edges
GEdge *firstEdge;
if (tempv.size() == 2){ // non periodic
firstEdge = (tempv.begin())->second;
@@ -256,3 +256,22 @@ SVector3 GEdgeCompound::firstDer(double par) const
getLocalParameter(par,iEdge,tLoc);
return _compound[iEdge]->firstDer(tLoc);
}
+
+void replaceMeshCompound(GFace *gf, std::list<GEdge*> &edges)
+{
+ std::list<GEdge*> e = gf->edges();
+ // replace edges by their compounds
+ std::set<GEdge*> mySet;
+ std::list<GEdge*>::iterator it = e.begin();
+ while(it != e.end()){
+ if((*it)->getCompound()){
+ mySet.insert((*it)->getCompound());
+ }
+ else{
+ mySet.insert(*it);
+ }
+ ++it;
+ }
+ edges.clear();
+ edges.insert(edges.begin(), mySet.begin(), mySet.end());
+}
diff --git a/Geo/GEdgeCompound.h b/Geo/GEdgeCompound.h
index a2c47fc78129880540f92704fd538e6b0c8cb91a..90bbe5fb7b9568e3e0c7e1e7f3f20c08b8e3f0e3 100644
--- a/Geo/GEdgeCompound.h
+++ b/Geo/GEdgeCompound.h
@@ -37,4 +37,6 @@ class GEdgeCompound : public GEdge {
std::vector<GEdge*> getCompounds() const { return _compound; }
};
+void replaceMeshCompound(GFace*, std::list<GEdge*>&);
+
#endif
diff --git a/Geo/GEdgeLoop.cpp b/Geo/GEdgeLoop.cpp
index 763cb786e9de6b756d61402f9949aefe77787580..728a312ea98be5285f9dc3b980f3c97fe9343a58 100644
--- a/Geo/GEdgeLoop.cpp
+++ b/Geo/GEdgeLoop.cpp
@@ -101,15 +101,13 @@ int GEdgeLoop::count(GEdge* ge) const
GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
{
-
-
- // Sometimes OCC puts a nasty degenerated edge in the middle of the wire ...
- // pushing it to front fixes the problem as it concerns gmsh
-
+ // Sometimes OCC puts a nasty degenerated edge in the middle of the
+ // wire ... pushing it to front fixes the problem as it concerns
+ // gmsh
std::list<GEdge*> wire;
std::vector<GEdge*> degenerated;
GEdge *degeneratedToInsert = 0;
- for (std::list<GEdge*>::const_iterator it = cwire.begin(); it != cwire.end() ; ++it){ // TEST
+ for (std::list<GEdge*>::const_iterator it = cwire.begin(); it != cwire.end(); ++it){
GEdge *ed = *it;
if (ed->degenerate(0))degenerated.push_back(ed);
else wire.push_back(ed);
@@ -129,25 +127,24 @@ GEdgeLoop::GEdgeLoop(const std::list<GEdge*> &cwire)
GEdgeSigned *prevOne = 0;
GEdgeSigned ges(0,0);
- // printf("go !!!\n");
while(wire.size()){
if (prevOne && degeneratedToInsert &&
degeneratedToInsert->getBeginVertex () == prevOne->getEndVertex()){
ges = GEdgeSigned(1,degeneratedToInsert);
degeneratedToInsert = 0;
- // printf("second degenerated edge inserted\n");
+ // printf("second degenerated edge inserted\n");
}
else ges = nextOne(prevOne, wire);
if(ges.getSign() == 0){ // oops
Msg::Error("Something wrong in edge loop of size=%d, no sign !", wire.size());
for (std::list<GEdge* >::iterator it = wire.begin(); it != wire.end(); it++){
- printf("GEdge=%d begin=%d end =%d \n", (*it)->tag(), (*it)->getBeginVertex()->tag(), (*it)->getEndVertex()->tag() );
+ Msg::Error("GEdge=%d begin=%d end =%d", (*it)->tag(),
+ (*it)->getBeginVertex()->tag(), (*it)->getEndVertex()->tag());
}
break;
}
prevOne = ⩾
- // ges.print();
+ // ges.print();
loop.push_back(ges);
}
}
-
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 1eb7f07e5f2c3db528331805d20b1e56ef7be4ff..4d33e5f1e2e526e9232aec0d52ce1e3885e3367c 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -51,8 +51,8 @@ static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssemb
}
}
-static void printBound(std::vector<MVertex*> &l, int tag){
- //print boundary
+static void printBound(std::vector<MVertex*> &l, int tag)
+{
char name[256];
sprintf(name, "myBOUND%d.pos", tag);
FILE * xyz = fopen(name,"w");
@@ -68,7 +68,6 @@ static void printBound(std::vector<MVertex*> &l, int tag){
static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
std::vector<double> &coord)
{
-
l.clear();
coord.clear();
@@ -130,7 +129,8 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
}
static void computeCGKernelPolygon(std::map<MVertex*,SPoint3> &coordinates,
- std::vector<MVertex*> &cavV, double &ucg, double &vcg)
+ std::vector<MVertex*> &cavV,
+ double &ucg, double &vcg)
{
ucg = 0.0;
vcg = 0.0;
@@ -237,13 +237,6 @@ static void myPolygon(std::vector<MElement*> &vTri, std::vector<MVertex*> &vPoly
}
}
- // printf("epoly.size=%d vTri=%d\n", ePoly.size(), vTri.size());
- // for(std::vector<MEdge>::iterator ite = ePoly.begin(); ite != ePoly.end(); ite++){
- // MVertex *vB = ite->getVertex(0);
- // MVertex *vE = ite->getVertex(1);
- // printf("VB=%d vE=%d \n", vB->getNum(), vE->getNum());
- // }
-
std::vector<MEdge>::iterator ite= ePoly.begin() ;
MVertex *vINIT = ite->getVertex(0);
vPoly.push_back(vINIT);
@@ -266,11 +259,6 @@ static void myPolygon(std::vector<MElement*> &vTri, std::vector<MVertex*> &vPoly
else ite++;
}
}
-
- // printf("epoly.size=%d vTri=%d, cavV.size =%d\n", ePoly.size(), vTri.size(), vPoly.size());
- // for(std::vector<MVertex*>::iterator itv = vPoly.begin(); itv != vPoly.end(); itv++){
- // printf("VV=%d \n", (*itv)->getNum());
- // }
}
bool checkCavity(std::vector<MElement*> &vTri, std::map<MVertex*, SPoint2> &vCoord)
@@ -351,27 +339,26 @@ void GFaceCompound::fillNeumannBCS() const
MVertex *v0 = orderedLoop[i];
MVertex *v1 = (i==nbLoop-1) ? orderedLoop[0]: orderedLoop[i+1];
- //loopfillTris.push_back(new MTriangle(v0,v1, c));
-
- // MVertex *v2 = new MVertex(.5*(v0->x()+c->x()), .5*(v0->y()+c->y()), .5*(v0->z()+c->z()));
- // MVertex *v3 = new MVertex(.5*(v1->x()+c->x()), .5*(v1->y()+c->y()), .5*(v1->z()+c->z()));
- // fillNodes.insert(c); fillNodes.insert(v2); fillNodes.insert(v3);
- // loopfillTris.push_back(new MTriangle(v0,v2,v3));
- // loopfillTris.push_back(new MTriangle(v2,c, v3));
- // loopfillTris.push_back(new MTriangle(v0,v3, v1));
-
- double k = 1/3.; double kk = 2/3.;
- MVertex *v2 = new MVertex(kk*v0->x()+k*c->x(), kk*v0->y()+k*c->y(),kk*v0->z()+k*c->z());
- MVertex *v3 = new MVertex(kk*v1->x()+k*c->x(), kk*v1->y()+k*c->y(),kk*v1->z()+k*c->z());
- MVertex *v4 = new MVertex(k*v0->x()+kk*c->x(), k*v0->y()+kk*c->y(),k*v0->z()+kk*c->z());
- MVertex *v5 = new MVertex(k*v1->x()+kk*c->x(), k*v1->y()+kk*c->y(),k*v1->z()+kk*c->z());
+ double k = 1. / 3.; double kk = 2. / 3.;
+ MVertex *v2 = new MVertex(kk * v0->x() + k * c->x(),
+ kk * v0->y() + k * c->y(),
+ kk * v0->z() + k * c->z());
+ MVertex *v3 = new MVertex(kk * v1->x() + k * c->x(),
+ kk * v1->y() + k * c->y(),
+ kk * v1->z() + k * c->z());
+ MVertex *v4 = new MVertex(k * v0->x() + kk * c->x(),
+ k * v0->y() + kk * c->y(),
+ k * v0->z() + kk * c->z());
+ MVertex *v5 = new MVertex(k * v1->x() + kk * c->x(),
+ k * v1->y() + kk * c->y(),
+ k * v1->z() + kk * c->z());
fillNodes.insert(c); fillNodes.insert(v2); fillNodes.insert(v3);
fillNodes.insert(v4); fillNodes.insert(v5);
- loopfillTris.push_back(new MTriangle(v0,v2,v3));
- loopfillTris.push_back(new MTriangle(v2,v5,v3));
- loopfillTris.push_back(new MTriangle(v2,v4,v5));
- loopfillTris.push_back(new MTriangle(v4,c,v5));
- loopfillTris.push_back(new MTriangle(v0,v3,v1));
+ loopfillTris.push_back(new MTriangle(v0, v2, v3));
+ loopfillTris.push_back(new MTriangle(v2, v5, v3));
+ loopfillTris.push_back(new MTriangle(v2, v4, v5));
+ loopfillTris.push_back(new MTriangle(v4, c, v5));
+ loopfillTris.push_back(new MTriangle(v0, v3, v1));
}
}
@@ -379,10 +366,12 @@ void GFaceCompound::fillNeumannBCS() const
//check normal orientations of loopfillTris
bool invertTris = false;
std::map<MEdge, std::set<MTriangle*>, Less_Edge > edge2tris;
- for(std::list<MTriangle*>::iterator t= loopfillTris.begin(); t!=loopfillTris.end(); t++){
+ for(std::list<MTriangle*>::iterator t = loopfillTris.begin();
+ t != loopfillTris.end(); t++){
for (int j = 0; j < 3; j++){
MEdge me = (*t)->getEdge(j);
- std::map<MEdge, std::set<MTriangle*, std::less<MTriangle*> >, Less_Edge >::iterator it = edge2tris.find(me);
+ std::map<MEdge, std::set<MTriangle*, std::less<MTriangle*> >,
+ Less_Edge >::iterator it = edge2tris.find(me);
if (it == edge2tris.end()) {
std::set<MTriangle*, std::less<MTriangle*> > mySet;
mySet.insert(*t);
@@ -402,7 +391,8 @@ void GFaceCompound::fillNeumannBCS() const
MTriangle *t = (*itf)->triangles[i];
for (int j = 0; j < 3; j++){
MEdge me = t->getEdge(j);
- std::map<MEdge, std::set<MTriangle*>, Less_Edge >::iterator it = edge2tris.find(me);
+ std::map<MEdge, std::set<MTriangle*>, Less_Edge >::iterator it =
+ edge2tris.find(me);
if(it != edge2tris.end()){
t->getEdgeInfo(me, iE,si);
MTriangle* t2 = *((it->second).begin());
@@ -416,12 +406,12 @@ void GFaceCompound::fillNeumannBCS() const
}
}
if (invertTris){
- for (std::list<MTriangle*>::iterator it = loopfillTris.begin(); it !=loopfillTris.end(); it++ )
+ for (std::list<MTriangle*>::iterator it = loopfillTris.begin();
+ it != loopfillTris.end(); it++ )
(*it)->revert();
}
fillTris.insert(fillTris.begin(),loopfillTris.begin(),loopfillTris.end());
-
}
//printing
@@ -432,7 +422,8 @@ void GFaceCompound::fillNeumannBCS() const
sprintf(name, "fillTris-%d.pos", (*itf)->tag());
FILE * ftri = fopen(name,"w");
fprintf(ftri,"View \"\"{\n");
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin();
+ it2 !=fillTris.end(); it2++ ){
MTriangle *t = (*it2);
fprintf(ftri,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
@@ -521,7 +512,6 @@ bool GFaceCompound::checkOverlap(std::vector<MVertex *> &vert) const
// the mapped triangles have the same normal orientation
bool GFaceCompound::checkOrientation(int iter) const
{
-
std::list<GFace*>::const_iterator it = _compound.begin();
double a_old = 0.0, a_new=0.0;
bool oriented = true;
@@ -611,7 +601,6 @@ void GFaceCompound::one2OneMap() const
bool GFaceCompound::parametrize() const
{
-
if (_compound.size() > 1) coherencePatches();
bool paramOK = true;
@@ -638,7 +627,8 @@ bool GFaceCompound::parametrize() const
// Multiscale Laplace parametrization
else if (_mapping == MULTISCALE){
std::vector<MElement*> _elements;
- for( std::list<GFace*>::const_iterator itt = _compound.begin(); itt != _compound.end(); ++itt)
+ for(std::list<GFace*>::const_iterator itt = _compound.begin();
+ itt != _compound.end(); ++itt)
for(unsigned int i = 0; i < (*itt)->triangles.size(); ++i)
_elements.push_back((*itt)->triangles[i]);
multiscaleLaplace multiLaplace(_elements, coordinates);
@@ -686,7 +676,7 @@ bool GFaceCompound::parametrize() const
buildOct();
if (_mapping != RBF){
- if (!checkOrientation(0) ){
+ if (!checkOrientation(0)){
Msg::Info("### Flipping: parametrization switched to convex combination map");
coordinates.clear();
Octree_Delete(oct);
@@ -709,8 +699,7 @@ bool GFaceCompound::parametrize() const
void GFaceCompound::getBoundingEdges()
{
-
- for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
+ for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
(*it)->setCompound(this);
}
@@ -751,7 +740,6 @@ void GFaceCompound::getBoundingEdges()
maxSize = size;
}
}
-
}
}
@@ -759,7 +747,8 @@ double GFaceCompound::getSizeH() const
{
SBoundingBox3d bb;
std::vector<SPoint3> vertices;
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ for(std::set<MVertex *>::iterator itv = allNodes.begin();
+ itv != allNodes.end() ; ++itv){
SPoint3 pt((*itv)->x(),(*itv)->y(), (*itv)->z());
vertices.push_back(pt);
bb +=pt;
@@ -838,7 +827,6 @@ SOrientedBoundingBox GFaceCompound::obb_boundEdges(const std::list<GEdge* > &eli
return res;
}
-
void GFaceCompound::getUniqueEdges(std::set<GEdge*> &_unique)
{
_unique.clear();
@@ -973,31 +961,30 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
return SPoint2(it->second.x(),it->second.y());
}
else{
-
double tGlob, tLoc;
double tL, tR;
int iEdge;
- //getParameter Point
+ // getParameter Point
v->getParameter(0,tGlob);
- //find compound Edge
+ // find compound Edge
GEdgeCompound *gec = dynamic_cast<GEdgeCompound*>(v->onWhat());
if(gec){
-
- //compute local parameter on Edge
+ // compute local parameter on Edge
gec->getLocalParameter(tGlob,iEdge,tLoc);
std::vector<GEdge*> gev = gec->getCompounds();
GEdge *ge = gev[iEdge];
- //left and right vertex of the Edge
+ // left and right vertex of the Edge
MVertex *v0 = ge->getBeginVertex()->mesh_vertices[0];
MVertex *v1 = ge->getEndVertex()->mesh_vertices[0];
std::map<MVertex*,SPoint3>::iterator itL = coordinates.find(v0);
std::map<MVertex*,SPoint3>::iterator itR = coordinates.find(v1);
- //for the Edge, find the left and right vertices of the initial 1D mesh and interpolate to find (u,v)
+ // for the Edge, find the left and right vertices of the initial
+ // 1D mesh and interpolate to find (u,v)
MVertex *vL = v0;
MVertex *vR = v1;
double tB = ge->parBounds(0).low();
@@ -1033,21 +1020,21 @@ SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
tR = tE;
}
- //Linear interpolation between tL and tR
+ // linear interpolation between tL and tR
double uloc, vloc;
uloc = itL->second.x() + (tLoc-tL)/(tR-tL) * (itR->second.x()-itL->second.x());
vloc = itL->second.y() + (tLoc-tL)/(tR-tL) * (itR->second.y()-itL->second.y());
return SPoint2(uloc,vloc);
}
}
+
// never here
return SPoint2(0, 0);
}
void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
{
-
-linearSystem<double> *_lsys = 0;
+ linearSystem<double> *_lsys = 0;
#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
_lsys = new linearSystemPETSc<double>;
#elif defined(HAVE_GMM) && !defined(HAVE_TAUCS)
@@ -1137,19 +1124,13 @@ linearSystem<double> *_lsys = 0;
}
delete _lsys;
-
}
bool GFaceCompound::parametrize_conformal_spectral() const
{
#if !defined(HAVE_PETSC) && !defined(HAVE_SLEPC)
-
- //Msg::Error("Gmsh should be compiled with petsc and slepc for using the conformal map.");
- //Msg::Error("Switch to harmonic map or see doc on the wiki for installing petsc and slepc:");
- //Msg::Error("https://geuz.org/trac/gmsh/wiki/STLRemeshing (username:gmsh,passwd:gmsh)");
Msg::Warning("Slepc not installed: parametrization switched to 'FE conformal' map");
- return parametrize_conformal(0,NULL,NULL);
-
+ return parametrize_conformal(0, NULL, NULL);
#else
linearSystem <double> *lsysA = new linearSystemPETSc<double>;
@@ -1228,7 +1209,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
}
// FIXME: force non-hermitian. For some reason (roundoff errors?)
- // on some machines et for some meshes slepc complains about bad IP
+ // on some machines and for some meshes slepc complains about bad IP
// norm otherwise
eigenSolver eig(&myAssembler, "B" , "A", false);
bool converged = eig.solve(1, "largest");
@@ -1259,22 +1240,20 @@ bool GFaceCompound::parametrize_conformal_spectral() const
return parametrize_conformal(0,NULL,NULL);
}
- std::vector<MVertex *> vert;
- bool hasOverlap = checkOverlap(vert);
- if (hasOverlap){
- Msg::Warning("!!! Overlap: parametrization switched to 'FE conformal' map");
- printStuff(3);
- return hasOverlap = parametrize_conformal(0, vert[0], vert[1]);
- }
-
- return hasOverlap;
-
+ std::vector<MVertex *> vert;
+ bool hasOverlap = checkOverlap(vert);
+ if (hasOverlap){
+ Msg::Warning("!!! Overlap: parametrization switched to 'FE conformal' map");
+ printStuff(3);
+ return hasOverlap = parametrize_conformal(0, vert[0], vert[1]);
+ }
+
+ return hasOverlap;
#endif
}
bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) const
{
-
linearSystem<double> *_lsys = 0;
#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
_lsys = new linearSystemPETSc<double>;
@@ -1290,8 +1269,8 @@ bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) co
dofManager<double> myAssembler(_lsys);
- if (!v1) v1 = _ordered[0];
- if (!v2) v2 = _ordered[(int)ceil((double)_ordered.size()/2.)];
+ if(!v1) v1 = _ordered[0];
+ if(!v2) v2 = _ordered[(int)ceil((double)_ordered.size()/2.)];
myAssembler.fixVertex(v1, 0, 1, 1.);
myAssembler.fixVertex(v1, 0, 2, 0.);
myAssembler.fixVertex(v2, 0, 1, -1.);
@@ -1309,7 +1288,8 @@ bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) co
myAssembler.numberVertex(t->getVertex(2), 0, 2);
}
}
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin();
+ it2 !=fillTris.end(); it2++){
MTriangle *t = (*it2);
myAssembler.numberVertex(t->getVertex(0), 0, 1);
myAssembler.numberVertex(t->getVertex(1), 0, 1);
@@ -1333,7 +1313,8 @@ bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) co
cross21.addToMatrix(myAssembler, &se);
}
}
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin();
+ it2 != fillTris.end(); it2++ ){
SElement se((*it2));
laplace1.addToMatrix(myAssembler, &se);
laplace2.addToMatrix(myAssembler, &se);
@@ -1345,7 +1326,8 @@ bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) co
_lsys->systemSolve();
Msg::Debug("System solved");
- for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ for(std::set<MVertex *>::iterator itv = allNodes.begin();
+ itv != allNodes.end() ; ++itv){
MVertex *v = *itv;
double value1, value2;
myAssembler.getDofValue(v, 0, 1, value1);
@@ -1355,19 +1337,19 @@ bool GFaceCompound::parametrize_conformal(int iter, MVertex *v1, MVertex *v2) co
delete _lsys;
-
- //check for overlap and compute new mapping with new pinned vertices
+ // check for overlap and compute new mapping with new pinned
+ // vertices
std::vector<MVertex *> vert;
bool hasOverlap = checkOverlap(vert);
if (hasOverlap && iter < 3){
- printf("**********Loop FE conformal iter (%d) v1=%d v2=%d \n", iter, vert[0]->getNum(), vert[1]->getNum());
- printStuff(100+iter);
- return hasOverlap = parametrize_conformal(iter+1, vert[0],vert[1]);
+ Msg::Info("Loop FE conformal iter (%d) v1=%d v2=%d", iter,
+ vert[0]->getNum(), vert[1]->getNum());
+ printStuff(100+iter);
+ return hasOverlap = parametrize_conformal(iter+1, vert[0],vert[1]);
}
else{
return hasOverlap;
}
-
}
void GFaceCompound::computeNormals(std::map<MVertex*,SVector3> &normals) const
@@ -1385,7 +1367,7 @@ void GFaceCompound::computeNormals() const
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
MTriangle *t = (*it)->triangles[i];
t->getJacobian(0, 0, 0, J);
- SVector3 n (J[2][0],J[2][1],J[2][2]);
+ SVector3 n (J[2][0], J[2][1], J[2][2]);
n.normalize();
for(int j = 0; j < 3; j++){
std::map<MVertex*, SVector3>::iterator itn = _normals.find(t->getVertex(j));
@@ -1400,25 +1382,22 @@ void GFaceCompound::computeNormals() const
double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
{
-
- if(!oct) parametrize();
- if(trivial()) {
+ if(!oct) parametrize();
+ if(trivial()) {
return (*(_compound.begin()))->curvatureMax(param);
- }
+ }
double U, V;
GFaceCompoundTriangle *lt;
getTriangle(param.x(), param.y(), <, U,V);
- if(!lt) {
- return 0.0;
- }
+ if(!lt) return 0.0;
- if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface) {
+ if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface) {
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
return lt->gf->curvatureMax(pv);
}
- else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
+ else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
Curvature& curvature = Curvature::getInstance();
if( !Curvature::valueAlreadyComputed() ) {
Msg::Info("Need to compute discrete curvature for isotropic remesh (in GFace)");
@@ -1433,33 +1412,30 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
return 0.;
}
+
double GFaceCompound::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
- double *curvMax, double *curvMin) const
+ double *curvMax, double *curvMin) const
{
-
- if(!oct) parametrize();
- if(trivial()) {
- return (*(_compound.begin()))->curvatures(param, dirMax,dirMin, curvMax, curvMin);
- }
+ if(!oct) parametrize();
+ if(trivial()) {
+ return (*(_compound.begin()))->curvatures(param, dirMax,dirMin, curvMax, curvMin);
+ }
double U, V;
GFaceCompoundTriangle *lt;
getTriangle(param.x(), param.y(), <, U,V);
- if(!lt) {
- return 0.0;
- }
+ if(!lt) return 0.0;
- if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface) {
+ if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface){
SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
return lt->gf->curvatures(pv, dirMax,dirMin, curvMax, curvMin);
}
-
- else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
+ else if (lt->gf->geomType() == GEntity::DiscreteSurface){
Curvature& curvature = Curvature::getInstance();
if( !Curvature::valueAlreadyComputed() ) {
Msg::Info("Need to compute discrete curvature for anisotropic remesh (in GFace)");
- Curvature::typeOfCurvature type = Curvature::RUSIN;//RBF
+ Curvature::typeOfCurvature type = Curvature::RUSIN; //RBF
curvature.computeCurvature(model(), type);
}
@@ -1470,28 +1446,25 @@ double GFaceCompound::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector
curvature.triangleNodalValuesAndDirections(lt->tri, dMax, dMin, cMax, cMin, 0);
//curvature.triangleNodalValuesAndDirections(lt->tri, dMax, dMin, cMax, cMin, 1);
- * dirMax = (1-U-V)*dMax[0] + U*dMax[1] + V*dMax[2];
- * dirMin = (1-U-V)*dMin[0] + U*dMin[1] + V*dMin[2];
- * curvMax = (1-U-V)*cMax[0] + U*cMax[1] + V*cMax[2];
- * curvMin = (1-U-V)*cMin[0] + U*cMin[1] + V*cMin[2];
+ *dirMax = (1-U-V)*dMax[0] + U*dMax[1] + V*dMax[2];
+ *dirMin = (1-U-V)*dMin[0] + U*dMin[1] + V*dMin[2];
+ *curvMax = (1-U-V)*cMax[0] + U*cMax[1] + V*cMax[2];
+ *curvMin = (1-U-V)*cMin[0] + U*cMin[1] + V*cMin[2];
- return * curvMax;
+ return *curvMax;
}
-
return 0.;
}
+
double GFaceCompound::locCurvature(MTriangle *t, double u, double v) const
{
-
SVector3 n1 = _normals[t->getVertex(0)];
SVector3 n2 = _normals[t->getVertex(1)];
SVector3 n3 = _normals[t->getVertex(2)];
double val[9] = {n1.x(), n2.x(), n3.x(),
n1.y(), n2.y(), n3.y(),
n1.z(), n2.z(), n3.z()};
-
return fabs(t->interpolateDiv(val, u, v, 3));
-
}
SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p, bool onSurface) const
@@ -1502,12 +1475,10 @@ SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p, bool onSurface) const
SPoint3 sp = it->second;
return SPoint2(sp.x(), sp.y());
-
}
GPoint GFaceCompound::point(double par1, double par2) const
{
-
if(trivial()){
return (*(_compound.begin()))->point(par1,par2);
}
@@ -1544,29 +1515,23 @@ GPoint GFaceCompound::point(double par1, double par2) const
const bool LINEARMESH = true; //false
if(LINEARMESH){
-
- //linear Lagrange mesh
- //-------------------------
+ // linear Lagrange mesh
p = lt->v1*(1.-U-V) + lt->v2*U + lt->v3*V;
return GPoint(p.x(),p.y(),p.z(),this,par);
-
}
else{
-
- //curved PN triangle
- //-------------------------
- printf("normals size=%d vertex=%d \n", (int)_normals.size(), lt->tri->getVertex(0)->getNum());
+ // curved PN triangle
const SVector3 n1 = _normals[lt->tri->getVertex(0)];
const SVector3 n2 = _normals[lt->tri->getVertex(1)];
const SVector3 n3 = _normals[lt->tri->getVertex(2)];
- SVector3 b300,b030,b003;
- SVector3 b210,b120,b021,b012,b102,b201,E,VV,b111;
+ SVector3 b300, b030, b003;
+ SVector3 b210, b120, b021, b012, b102, b201, E, VV, b111;
b300 = lt->v1;
b030 = lt->v2;
b003 = lt->v3;
-
- //tagged PN triangles (sigma=1)
+
+ // tagged PN triangles (sigma=1)
double theta = 0.0;
SVector3 d1 = lt->v1+.33*(1-theta)*(lt->v2-lt->v1);
SVector3 d2 = lt->v2+.33*(1-theta)*(lt->v1-lt->v2);
@@ -1604,7 +1569,6 @@ GPoint GFaceCompound::point(double par1, double par2) const
Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
{
-
if(!oct) parametrize();
if(trivial())
@@ -1620,27 +1584,27 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
SVector3 dXdu, dXdv ;
bool conv = _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
return Pair<SVector3, SVector3>(dXdu,dXdv);
- }
+ }
double mat[2][2] = {{lt->p2.x() - lt->p1.x(), lt->p3.x() - lt->p1.x()},
{lt->p2.y() - lt->p1.y(), lt->p3.y() - lt->p1.y()}};
double inv[2][2];
double det = inv2x2(mat,inv);
if (!det && _mapping == RBF){
- double x, y, z;
- SVector3 dXdu, dXdv ;
- bool conv = _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
- return Pair<SVector3, SVector3>(dXdu,dXdv);
+ double x, y, z;
+ SVector3 dXdu, dXdv ;
+ bool conv = _rbf->UVStoXYZ(param.x(), param.y(), x,y,z, dXdu, dXdv);
+ return Pair<SVector3, SVector3>(dXdu,dXdv);
}
-
+
SVector3 dXdxi(lt->v2 - lt->v1);
SVector3 dXdeta(lt->v3 - lt->v1);
-
+
SVector3 dXdu(dXdxi * inv[0][0] + dXdeta * inv[1][0]);
SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]);
-
- return Pair<SVector3, SVector3>(dXdu,dXdv);
-}
+
+ return Pair<SVector3, SVector3>(dXdu, dXdv);
+}
void GFaceCompound::secondDer(const SPoint2 ¶m,
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const
@@ -1648,7 +1612,6 @@ void GFaceCompound::secondDer(const SPoint2 ¶m,
if(!oct) parametrize();
//leave debug here (since outputScalarField calls curvatureDiv)
Msg::Debug("Computation of the second derivatives is not implemented for compound faces");
-
}
static void GFaceCompoundBB(void *a, double*mmin, double*mmax)
@@ -1681,7 +1644,7 @@ static void GFaceCompoundCentroid(void *a, double*c)
static int GFaceCompoundInEle(void *a, double*c)
{
GFaceCompoundTriangle *t = (GFaceCompoundTriangle *)a;
- double M[2][2],R[2],X[2];
+ double M[2][2], R[2], X[2];
const double eps = 1.e-8;
const SPoint3 p0 = t->p1;
const SPoint3 p1 = t->p2;
@@ -1692,7 +1655,7 @@ static int GFaceCompoundInEle(void *a, double*c)
M[1][1] = p2.y() - p0.y();
R[0] = (c[0] - p0.x());
R[1] = (c[1] - p0.y());
- sys2x2(M,R,X);
+ sys2x2(M, R, X);
if(X[0] > -eps && X[1] > -eps && 1. - X[0] - X[1] > -eps){
return 1;
}
@@ -1742,15 +1705,13 @@ void GFaceCompound::getTriangle(double u, double v,
M[1][1] = p2.y() - p0.y();
R[0] = (u - p0.x());
R[1] = (v - p0.y());
- sys2x2(M,R,X);
+ sys2x2(M, R, X);
_u = X[0];
_v = X[1];
}
-
void GFaceCompound::buildOct() const
{
-
SBoundingBox3d bb;
int count = 0;
std::list<GFace*>::const_iterator it = _compound.begin();
@@ -1793,16 +1754,22 @@ void GFaceCompound::buildOct() const
if((*it)->geomType() != GEntity::DiscreteSurface){
// take care of the seam !!!!
if (t->getVertex(0)->onWhat()->dim() == 2){
- reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(1),*it, _gfct[count].gfp1, _gfct[count].gfp2);
- reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(2),*it, _gfct[count].gfp1, _gfct[count].gfp3);
+ reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(1), *it,
+ _gfct[count].gfp1, _gfct[count].gfp2);
+ reparamMeshEdgeOnFace(t->getVertex(0), t->getVertex(2), *it,
+ _gfct[count].gfp1, _gfct[count].gfp3);
}
else if (t->getVertex(1)->onWhat()->dim() == 2){
- reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(0),*it, _gfct[count].gfp2, _gfct[count].gfp1);
- reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(2),*it, _gfct[count].gfp2, _gfct[count].gfp3);
+ reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(0), *it,
+ _gfct[count].gfp2, _gfct[count].gfp1);
+ reparamMeshEdgeOnFace(t->getVertex(1), t->getVertex(2), *it,
+ _gfct[count].gfp2, _gfct[count].gfp3);
}
else if (t->getVertex(2)->onWhat()->dim() == 2){
- reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(0),*it, _gfct[count].gfp3, _gfct[count].gfp1);
- reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(1),*it, _gfct[count].gfp3, _gfct[count].gfp2);
+ reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(0), *it,
+ _gfct[count].gfp3, _gfct[count].gfp1);
+ reparamMeshEdgeOnFace(t->getVertex(2), t->getVertex(1), *it,
+ _gfct[count].gfp3, _gfct[count].gfp2);
}
else {
reparamMeshVertexOnFace(t->getVertex(0), *it, _gfct[count].gfp1);
@@ -1811,27 +1778,25 @@ void GFaceCompound::buildOct() const
}
}
_gfct[count].v1 = SPoint3(t->getVertex(0)->x(), t->getVertex(0)->y(),
- t->getVertex(0)->z());
+ t->getVertex(0)->z());
_gfct[count].v2 = SPoint3(t->getVertex(1)->x(), t->getVertex(1)->y(),
- t->getVertex(1)->z());
+ t->getVertex(1)->z());
_gfct[count].v3 = SPoint3(t->getVertex(2)->x(), t->getVertex(2)->y(),
- t->getVertex(2)->z());
- _gfct[count].gf = *it;
- _gfct[count].tri = t;
+ t->getVertex(2)->z());
+ _gfct[count].gf = *it;
+ _gfct[count].tri = t;
Octree_Insert(&_gfct[count], oct);
- count ++;
+ count++;
}
}
nbT = count;
Octree_Arrange(oct);
printStuff();
-
}
bool GFaceCompound::checkTopology() const
{
-
if (_mapping == RBF) return true;
//fixme tristan
@@ -1856,12 +1821,13 @@ bool GFaceCompound::checkTopology() const
Msg::Warning("Wrong topology: Genus=%d, Nb boundaries=%d, AR=%g", G, Nb, H/D);
if (_allowPartition){
Msg::Info("-----------------------------------------------------------");
- Msg::Info("--- Split surface %d in %d parts with Multilevel Mesh partitioner", tag(), nbSplit);
+ Msg::Info("--- Split surface %d in %d parts with Multilevel Mesh partitioner",
+ tag(), nbSplit);
}
else{
- Msg::Error("For remeshing your geometry, you should enable the automatic remeshing algorithm.");
- Msg::Error("Add 'Mesh.RemeshAlgorithm=1;' in your geo file or through the Fltk window (Options > Mesh > General)");
- Msg::Exit(0);
+ Msg::Fatal("For remeshing your geometry, you should enable the automatic "
+ "remeshing algorithm. Add 'Mesh.RemeshAlgorithm=1;' in your "
+ "geo file or through the Fltk window (Options > Mesh > General)");
}
}
else if (G == 0 && AR > AR_MAX){
@@ -1870,25 +1836,28 @@ bool GFaceCompound::checkTopology() const
if (_allowPartition == 1){
nbSplit = -2;
Msg::Info("-----------------------------------------------------------");
- Msg::Info("--- Split surface %d in 2 parts with Laplacian Mesh partitioner", tag());
+ Msg::Info("--- Split surface %d in 2 parts with Laplacian Mesh partitioner",
+ tag());
}
else if (_allowPartition == 2){
nbSplit = 2;
Msg::Info("-----------------------------------------------------------");
- Msg::Info("--- Split surface %d in %d parts with Multilevel Mesh partitioner", tag(), nbSplit);
+ Msg::Info("--- Split surface %d in %d parts with Multilevel Mesh partitioner",
+ tag(), nbSplit);
}
else if (_allowPartition == 0){
- Msg::Warning("The geometrical aspect ratio of your geometry is quite high.");
- Msg::Warning("You should enable partitioning of the mesh by activating the automatic remeshin algorithm.");
- Msg::Warning("Add 'Mesh.RemeshAlgorithm=1;' in your geo file or through the Fltk window (Options > Mesh > General)");
+ Msg::Warning("The geometrical aspect ratio of your geometry is quite high. "
+ "You should enable partitioning of the mesh by activating the "
+ "automatic remeshin algorithm. Add 'Mesh.RemeshAlgorithm=1;' "
+ "in your geo file or through the Fltk window (Options > Mesh > "
+ "General)");
}
}
else{
Msg::Debug("Correct topology: Genus=%d and Nb boundaries=%d, AR=%g", G, Nb, H/D);
}
-
+
return correctTopo;
-
}
double GFaceCompound::checkAspectRatio() const
@@ -1949,14 +1918,13 @@ double GFaceCompound::checkAspectRatio() const
}
return AR;
-
}
void GFaceCompound::coherencePatches() const
{
-
if (_mapping == RBF) return;
- Msg::Info("Re-orient all %d compound patches normals coherently", _compound.size());
+ Msg::Info("Re-orient all %d compound patches normals coherently",
+ _compound.size());
std::map<MEdge, std::set<MElement*>, Less_Edge > edge2elems;
std::vector<MElement*> allElems;
@@ -1967,7 +1935,8 @@ void GFaceCompound::coherencePatches() const
allElems.push_back(t);
for (int j = 0; j < t->getNumEdges(); j++){
MEdge me = t->getEdge(j);
- std::map<MEdge, std::set<MElement*, std::less<MElement*> >, Less_Edge >::iterator it = edge2elems.find(me);
+ std::map<MEdge, std::set<MElement*, std::less<MElement*> >,
+ Less_Edge >::iterator it = edge2elems.find(me);
if (it == edge2elems.end()) {
std::set<MElement*, std::less<MElement*> > mySet;
mySet.insert(t);
@@ -1993,9 +1962,11 @@ void GFaceCompound::coherencePatches() const
for (int j = 0; j < t->getNumEdges(); j++){
MEdge me = t->getEdge(j);
t->getEdgeInfo(me, iE,si);
- std::map<MEdge, std::set<MElement*>, Less_Edge >::iterator it = edge2elems.find(me);
+ std::map<MEdge, std::set<MElement*>, Less_Edge >::iterator it =
+ edge2elems.find(me);
std::set<MElement*, std::less<MElement*> > mySet = it->second;
- for(std::set<MElement*, std::less<MElement*> >::iterator itt = mySet.begin(); itt != mySet.end(); itt++){
+ for(std::set<MElement*, std::less<MElement*> >::iterator itt = mySet.begin();
+ itt != mySet.end(); itt++){
if (*itt != t){
(*itt)->getEdgeInfo(me,iE2,si2);
if(si == si2) (*itt)->revert();
@@ -2017,7 +1988,8 @@ void GFaceCompound::coherenceNormals()
MElement *t = getMeshElement(i);
for (int j = 0; j < t->getNumEdges(); j++){
MEdge me = t->getEdge(j);
- std::map<MEdge, std::set<MElement*, std::less<MElement*> >, Less_Edge >::iterator it = edge2elems.find(me);
+ std::map<MEdge, std::set<MElement*, std::less<MElement*> >,
+ Less_Edge >::iterator it = edge2elems.find(me);
if (it == edge2elems.end()) {
std::set<MElement*, std::less<MElement*> > mySet;
mySet.insert(t);
@@ -2042,9 +2014,11 @@ void GFaceCompound::coherenceNormals()
for (int j = 0; j < t->getNumEdges(); j++){
MEdge me = t->getEdge(j);
t->getEdgeInfo(me, iE,si);
- std::map<MEdge, std::set<MElement*>, Less_Edge >::iterator it = edge2elems.find(me);
+ std::map<MEdge, std::set<MElement*>, Less_Edge >::iterator it =
+ edge2elems.find(me);
std::set<MElement*, std::less<MElement*> > mySet = it->second;
- for(std::set<MElement*, std::less<MElement*> >::iterator itt = mySet.begin(); itt != mySet.end(); itt++){
+ for(std::set<MElement*, std::less<MElement*> >::iterator itt = mySet.begin();
+ itt != mySet.end(); itt++){
if (*itt != t){
(*itt)->getEdgeInfo(me,iE2,si2);
if(si == si2) (*itt)->revert();
@@ -2069,7 +2043,6 @@ void GFaceCompound::buildAllNodes() const
}
}
}
-
}
int GFaceCompound::genusGeom() const
@@ -2123,15 +2096,14 @@ void GFaceCompound::printStuff(int iNewton) const
//FILE * xyzc = fopen(name6,"w");
//FILE * uvm = fopen(name7,"w");
-
- //fprintf(uva,"View \"\"{\n");
- fprintf(uvx,"View \"\"{\n");
- fprintf(uvy,"View \"\"{\n");
- fprintf(uvz,"View \"\"{\n");
- fprintf(xyzu,"View \"\"{\n");
- fprintf(xyzv,"View \"\"{\n");
- //fprintf(xyzc,"View \"\"{\n");
- //fprintf(uvm,"View \"\"{\n");
+ //fprintf(uva, "View \"\"{\n");
+ fprintf(uvx, "View \"\"{\n");
+ fprintf(uvy, "View \"\"{\n");
+ fprintf(uvz, "View \"\"{\n");
+ fprintf(xyzu, "View \"\"{\n");
+ fprintf(xyzv, "View \"\"{\n");
+ //fprintf(xyzc, "View \"\"{\n");
+ //fprintf(uvm, "View \"\"{\n");
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -2142,12 +2114,13 @@ void GFaceCompound::printStuff(int iNewton) const
coordinates.find(t->getVertex(1));
std::map<MVertex*,SPoint3>::const_iterator it2 =
coordinates.find(t->getVertex(2));
- fprintf(xyzv,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(xyzv, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
it0->second.y(),it1->second.y(),it2->second.y());
- fprintf(xyzu,"ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
+ fprintf(xyzu, "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,"
+ "%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
@@ -2171,17 +2144,23 @@ void GFaceCompound::printStuff(int iNewton) const
double a_2D = fabs(triangle_area(q0, q1, q2));
double area = (a_3D/a_2D); //*(a_3D/a_2D);
- // Pair<SVector3, SVector3> der = this->firstDer(SPoint2( it0->second.x(), it0->second.y()));
+ // Pair<SVector3, SVector3> der = this->firstDer(SPoint2(it0->second.x(),
+ // it0->second.y()));
// double metric0e = dot(der.first(), der.first());
- // double metric0f = dot(der.second()*(1./norm(der.second())), der.first()*(1./norm(der.first())));
+ // double metric0f = dot(der.second()*(1./norm(der.second())),
+ // der.first()*(1./norm(der.first())));
// double metric0g = dot(der.second(), der.second());
- // Pair<SVector3, SVector3> der1 = this->firstDer(SPoint2( it1->second.x(), it1->second.y()));
+ // Pair<SVector3, SVector3> der1 = this->firstDer(SPoint2(it1->second.x(),
+ // it1->second.y()));
// double metric1e = dot(der1.first(), der1.first());
- // double metric1f = dot(der1.second()*(1/norm(der1.second())), der1.first()*(1./norm(der1.first())));
+ // double metric1f = dot(der1.second()*(1/norm(der1.second())),
+ // der1.first()*(1./norm(der1.first())));
// double metric1g = dot(der1.second(), der1.second());
- // Pair<SVector3, SVector3> der2 = this->firstDer(SPoint2( it2->second.x(), it2->second.y()));
+ // Pair<SVector3, SVector3> der2 = this->firstDer(SPoint2(it2->second.x(),
+ // it2->second.y()));
// double metric2e = dot(der2.first(), der2.first());
- // double metric2f = dot(der2.second()*(1./norm(der2.second())), der2.first()*(1./norm(der2.first())));
+ // double metric2f = dot(der2.second()*(1./norm(der2.second())),
+ // der2.first()*(1./norm(der2.first())));
// double metric2g = dot(der2.second(), der2.second());
// double mat0[2][2], eig0[2];
@@ -2201,29 +2180,30 @@ void GFaceCompound::printStuff(int iNewton) const
// double disp1 = sqrt(.5*(eig1[0]*eig1[0]+ (eig1[1]*eig1[1])));
// double disp2 = sqrt(.5*(eig2[0]*eig2[0]+ (eig2[1]*eig2[1])));
// double mdisp = .333*(disp0+disp1+disp2);
- // fprintf(uva,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ // fprintf(uva, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
// it0->second.x(), it0->second.y(), 0.0,
// it1->second.x(), it1->second.y(), 0.0,
// it2->second.x(), it2->second.y(), 0.0,
// area, area, area);
- // fprintf(uvm,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ // fprintf(uvm, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
// it0->second.x(), it0->second.y(), 0.0,
// it1->second.x(), it1->second.y(), 0.0,
// it2->second.x(), it2->second.y(), 0.0,
// mdisp, mdisp, mdisp);
- fprintf(uvx,"ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
+ fprintf(uvx, "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,"
+ "%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
it0->second.x(), it0->second.y(), 0.0,
it1->second.x(), it1->second.y(), 0.0,
it2->second.x(), it2->second.y(), 0.0,
t->getVertex(0)->x(), t->getVertex(1)->x(), t->getVertex(2)->x());
- fprintf(uvy,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(uvy, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(), it0->second.y(), 0.0,
it1->second.x(), it1->second.y(), 0.0,
it2->second.x(), it2->second.y(), 0.0,
t->getVertex(0)->y(), t->getVertex(1)->y(), t->getVertex(2)->y());
- fprintf(uvz,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(uvz, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(), it0->second.y(), 0.0,
it1->second.x(), it1->second.y(), 0.0,
it2->second.x(), it2->second.y(), 0.0,
@@ -2250,21 +2230,3 @@ void GFaceCompound::printStuff(int iNewton) const
#endif
-void replaceMeshCompound(GFace *gf, std::list<GEdge*> &edges)
-{
- std::list<GEdge*> e = gf->edges();
- //Replace edges by their compounds
- std::set<GEdge*> mySet;
- std::list<GEdge*>::iterator it = e.begin();
- while(it != e.end()){
- if((*it)->getCompound()){
- mySet.insert((*it)->getCompound());
- }
- else{
- mySet.insert(*it);
- }
- ++it;
- }
- edges.clear();
- edges.insert(edges.begin(), mySet.begin(), mySet.end());
-}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index bd316637053f94a814f690fe1490d98308a47caa..74aef07aad073da7433d332211c2d9df4c72af7a 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -169,6 +169,4 @@ class GFaceCompound : public GFace {
#endif
-void replaceMeshCompound(GFace*,std::list<GEdge*>&);
-
#endif
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index 52f2521bcb1cf6b02ea6692900668e11776b101c..df203b3332322d0970a651fafc4d8a12374980a5 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -4,7 +4,6 @@
// bugs and problems to <gmsh@geuz.org>.
#include "GModelFactory.h"
-
#include "ListUtils.h"
#include "Context.h"
#include "GVertex.h"
@@ -15,7 +14,6 @@
#include "MLine.h"
#include "GModel.h"
-
GVertex *GeoFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
{
int num = gm->getMaxElementaryNumber(0) + 1;
@@ -79,9 +77,11 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
List_T *iListl = List_Create(nl, nl, sizeof(int));
if (nl > 1) {
const GEdge &e1 = *edges[i][1];
- if (e0.getEndVertex() != e1.getBeginVertex() && e0.getEndVertex() != e1.getEndVertex()) {
+ if (e0.getEndVertex() != e1.getBeginVertex() &&
+ e0.getEndVertex() != e1.getEndVertex()) {
frontV = e0.getEndVertex();
- if (e0.getBeginVertex() != e1.getBeginVertex() && e0.getBeginVertex() != e1.getEndVertex()) {
+ if (e0.getBeginVertex() != e1.getBeginVertex() &&
+ e0.getBeginVertex() != e1.getEndVertex()) {
Msg::Warning("Edges 0 and 1 not consecutive in line loop %d", i);
//return NULL;
}
@@ -183,7 +183,6 @@ GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > f
return gr;
}
-//---------------------------------------------------------------------------------
#if defined(HAVE_OCC)
#include "OCCIncludes.h"
#include "GModelIO_OCC.h"
@@ -215,7 +214,6 @@ GRegion* GeoFactory::addVolume (GModel *gm, std::vector<std::vector<GFace *> > f
#include <TColStd_HArray1OfReal.hxx>
#include <TColStd_HArray1OfInteger.hxx>
-
GVertex *OCCFactory::addVertex(GModel *gm, double x, double y, double z, double lc)
{
if (!gm->_occ_internals)
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index d6094e6dbcfa16b1b25aa109a746b3faacde1dee..e7aa84ef0e047b9424854c5e800909d18e6f6bc0 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -1,5 +1,11 @@
-#include "GRbf.h"
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by FIXME
+#include "GRbf.h"
#include <math.h>
#include <vector>
#include "fullMatrix.h"
@@ -15,7 +21,8 @@
#include <ANN/ANN.h>
#endif
-static void SphereBB(void *a, double*mmin, double*mmax){
+static void SphereBB(void *a, double*mmin, double*mmax)
+{
Sphere *t = (Sphere *)a;
mmin[0] = t->center.x()-t->radius;
mmin[1] = t->center.y()-t->radius;
@@ -25,14 +32,16 @@ static void SphereBB(void *a, double*mmin, double*mmax){
mmax[2] = t->center.z()+t->radius;
}
-static void SphereCentroid(void *a, double*c){
+static void SphereCentroid(void *a, double*c)
+{
Sphere *t = (Sphere *)a;
c[0] = t->center.x();
c[1] = t->center.y();
c[2] = t->center.z();
}
-static int SphereInEle(void *a, double*c){
+static int SphereInEle(void *a, double*c)
+{
Sphere *t = (Sphere *)a;
double dist = sqrt((c[0]-t->center.x())*(c[0]-t->center.x())+
(c[1]-t->center.y())*(c[1]-t->center.y())+
@@ -42,7 +51,9 @@ static int SphereInEle(void *a, double*c){
}
return 0;
}
-static void printNodes(std::set<MVertex *> myNodes){
+
+static void printNodes(std::set<MVertex *> myNodes)
+{
FILE * xyz = fopen("myNodes.pos","w");
fprintf(xyz,"View \"\"{\n");
for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end() ; ++itv){
@@ -53,9 +64,8 @@ static void printNodes(std::set<MVertex *> myNodes){
fclose(xyz);
}
-
-static void exportParametrizedMesh(fullMatrix<double> &UV, int nbNodes){
-
+static void exportParametrizedMesh(fullMatrix<double> &UV, int nbNodes)
+{
FILE *f = fopen ("UV.pos", "w");
fprintf(f,"View \" uv \" {\n");
@@ -68,27 +78,28 @@ static void exportParametrizedMesh(fullMatrix<double> &UV, int nbNodes){
fclose(f);
}
-GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVector3> _normals,
- std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes, bool _isLocal)
- : sBox(sizeBox), variableShapeParam(variableEps), radialFunctionIndex (rbfFun), _inUV(0), isLocal(_isLocal)
+GRbf::GRbf(double sizeBox, int variableEps, int rbfFun,
+ std::map<MVertex*, SVector3> _normals,
+ std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes, bool _isLocal)
+ : sBox(sizeBox), variableShapeParam(variableEps), radialFunctionIndex (rbfFun),
+ _inUV(0), isLocal(_isLocal)
{
-
#if defined (HAVE_ANN)
- XYZkdtree=0;
+ XYZkdtree = 0;
#endif
allCenters.resize(allNodes.size(),3);
double tol = 4.e-2*sBox;
if (isLocal) tol = 1.e-15;
- //remove duplicate vertices
- //add bc nodes
+ // remove duplicate vertices
+ // add bc nodes
for(unsigned int i = 0; i < bcNodes.size(); i++){
myNodes.insert(bcNodes[i]);
//if (bcNodes.size() > 20) i+=2;
}
- //then create Mvertex position
+ // then create Mvertex position
std::vector<MVertex*> vertices( allNodes.begin(), allNodes.end() );
MVertexPositionSet pos(vertices);
for(unsigned int i = 0; i < vertices.size(); i++){
@@ -99,18 +110,18 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
allCenters(i,2) = v->z()/sBox;
_mapAllV.insert(std::make_pair(v, i));
}
- //keep only no duplicate vertices
+ // keep only no duplicate vertices
for(unsigned int i = 0; i < vertices.size(); i++)
if(vertices[i]->getIndex()) myNodes.insert(vertices[i]);
- //initialize with points
+ // initialize with points
nbNodes= myNodes.size();
centers.resize(nbNodes,3);
normals.resize(nbNodes,3);
int index = 0;
double dist_min = 1.e6;
double dist_max = 1.e-6;
- for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end() ; ++itv){
+ for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end(); ++itv){
MVertex *v1 = *itv;
centers(index,0) = v1->x()/sBox;
centers(index,1) = v1->y()/sBox;
@@ -122,10 +133,11 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
normals(index,2) = itn->second.z();
}
_mapV.insert(std::make_pair(v1, index));
- for(std::set<MVertex *>::iterator itv2 = myNodes.begin(); itv2 !=myNodes.end() ; itv2++){
+ for(std::set<MVertex *>::iterator itv2 = myNodes.begin(); itv2 !=myNodes.end(); itv2++){
MVertex *v2 = *itv2;
- double dist = sqrt((v1->x()-v2->x())*(v1->x()-v2->x())+(v1->y()-v2->y())*(v1->y()-v2->y())
- +(v1->z()-v2->z())*(v1->z()-v2->z()))/sBox;
+ double dist = sqrt((v1->x()-v2->x())*(v1->x()-v2->x())+
+ (v1->y()-v2->y())*(v1->y()-v2->y())+
+ (v1->z()-v2->z())*(v1->z()-v2->z()))/sBox;
if (dist<dist_min && *itv != *itv2 && dist > 1.e-5) dist_min = dist;
}
index++;
@@ -148,10 +160,10 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
}
extendedX.resize(3*nbNodes,3);
-
}
-GRbf::~GRbf(){
+GRbf::~GRbf()
+{
#if defined (HAVE_ANN)
ANNpointArray XYZNodes = XYZkdtree->thePoints();
ANNpointArray UVNodes = UVkdtree->thePoints();
@@ -160,11 +172,10 @@ GRbf::~GRbf(){
delete XYZkdtree;
delete UVkdtree;
#endif
-
}
-void GRbf::buildXYZkdtree(){
-
+void GRbf::buildXYZkdtree()
+{
#if defined (HAVE_ANN)
ANNpointArray XYZnodes = annAllocPts(nbNodes, 3);
for(int i = 0; i < nbNodes; i++){
@@ -174,11 +185,10 @@ void GRbf::buildXYZkdtree(){
}
XYZkdtree = new ANNkd_tree(XYZnodes, nbNodes, 3);
#endif
-
}
-void GRbf::buildOctree(double radius){
-
+void GRbf::buildOctree(double radius)
+{
printf("building octree radius = %g \n", radius);
SBoundingBox3d bb;
for (int i= 0; i< nbNodes; i++)
@@ -215,15 +225,13 @@ void GRbf::buildOctree(double radius){
}
Octree_Delete(oct);
-
buildXYZkdtree();
-
}
-//compute curvature from level set
+// compute curvature from level set
void GRbf::curvatureRBF(const fullMatrix<double> &cntrs,
- fullMatrix<double> &curvature){
-
+ fullMatrix<double> &curvature)
+{
fullMatrix<double> extX, surf, sx,sy,sz, sxx,syy,szz, sxy,sxz,syz,sLap;
setup_level_set(cntrs,normals,extX, surf);
@@ -240,8 +248,8 @@ void GRbf::curvatureRBF(const fullMatrix<double> &cntrs,
}
void GRbf::computeCurvature(const fullMatrix<double> &cntrs,
- std::map<MVertex*, double> &rbf_curv){
-
+ std::map<MVertex*, double> &rbf_curv)
+{
fullMatrix<double> curvature(cntrs.size1(), 1);
curvatureRBF(cntrs, curvature);
@@ -259,8 +267,8 @@ void GRbf::computeCurvature(const fullMatrix<double> &cntrs,
}
void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
- std::map<MVertex*, double> &rbf_curv) {
-
+ std::map<MVertex*, double> &rbf_curv)
+{
int numNodes = cntrs.size1();
if(nodesInSphere.size() == 0) buildOctree(radius);
@@ -291,8 +299,8 @@ void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
}
-double GRbf::evalRadialFnDer (int p, double dx, double dy, double dz, double ep){
-
+double GRbf::evalRadialFnDer (int p, double dx, double dy, double dz, double ep)
+{
double r2 = dx*dx+dy*dy+dz*dz; //r^2
switch (radialFunctionIndex) {
case 0 : // GA
@@ -328,8 +336,8 @@ double GRbf::evalRadialFnDer (int p, double dx, double dy, double dz, double ep)
fullMatrix<double> GRbf::generateRbfMat(int p,
const fullMatrix<double> &nodes1,
- const fullMatrix<double> &nodes2) {
-
+ const fullMatrix<double> &nodes2)
+{
int m = nodes2.size1();
int n = nodes1.size1();
fullMatrix<double> rbfMat(m,n);
@@ -346,15 +354,14 @@ fullMatrix<double> GRbf::generateRbfMat(int p,
}
return rbfMat;
-
}
//DL matrix (B*inv A)
void GRbf::RbfOp(int p,
const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes,
- fullMatrix<double> &D) {
-
+ fullMatrix<double> &D)
+{
fullMatrix<double> rbfInvA, rbfMatB;
D.resize(nodes.size1(), cntrs.size1());
@@ -383,26 +390,26 @@ void GRbf::evalRbfDer(int p,
const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes,
const fullMatrix<double> &fValues,
- fullMatrix<double> &fApprox) {
-
+ fullMatrix<double> &fApprox)
+{
fApprox.resize(nodes.size1(),fValues.size2());
fullMatrix<double> D;
RbfOp(p,cntrs,nodes,D);
fApprox.gemm(D,fValues, 1.0, 0.0);
-
}
void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
fullMatrix<double> &level_set_nodes,
- fullMatrix<double> &level_set_funvals){
-
+ fullMatrix<double> &level_set_funvals)
+{
int numNodes = cntrs.size1();
int nTot = 3*numNodes;
double normFactor;
level_set_nodes.resize(nTot,3);
level_set_funvals.resize(nTot,1);
- fullMatrix<double> ONES(numNodes+1,1),sx(numNodes,1),sy(numNodes,1),sz(numNodes,1),norms(numNodes,3), cntrsPlus(numNodes+1,3);
+ fullMatrix<double> ONES(numNodes+1,1), sx(numNodes,1), sy(numNodes,1);
+ fullMatrix<double> sz(numNodes,1),norms(numNodes,3), cntrsPlus(numNodes+1,3);
//Computes the normal vectors to the surface at each node
for (int i=0;i<numNodes ; ++i){
@@ -446,13 +453,12 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
matAInv_nn.resize(nTot, nTot);
matAInv_nn = generateRbfMat(0,level_set_nodes,level_set_nodes);
matAInv_nn.invertInPlace();
-
}
-
void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
- fullMatrix<double> &Oper) {
+ fullMatrix<double> &Oper)
+{
isLocal = true;
int numNodes = cntrs.size1();
Oper.resize(numNodes,numNodes);
@@ -485,8 +491,8 @@ void GRbf::RbfLapSurface_local_projection(const fullMatrix<double> &cntrs,
void GRbf::RbfLapSurface_global_projection( const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
- fullMatrix<double> &Oper){
-
+ fullMatrix<double> &Oper)
+{
int numNodes = cntrs.size1();
Oper.resize(numNodes,numNodes);
@@ -536,8 +542,8 @@ void GRbf::RbfLapSurface_global_projection( const fullMatrix<double> &cntrs,
void GRbf::RbfLapSurface_local_CPM(bool isLow,
const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
- fullMatrix<double> &Oper){
-
+ fullMatrix<double> &Oper)
+{
isLocal = true;
int numNodes = cntrs.size1();
Oper.resize(3*numNodes,3*numNodes);
@@ -582,10 +588,11 @@ void GRbf::RbfLapSurface_local_CPM(bool isLow,
}
}
-//NEW METHOD #1 CPM GLOBAL HIGH
+// NEW METHOD #1 CPM GLOBAL HIGH
void GRbf::RbfLapSurface_global_CPM_high_2(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
- fullMatrix<double> &Oper){
+ fullMatrix<double> &Oper)
+{
Msg::Info("*** RBF ... building Laplacian operator");
int numNodes = cntrs.size1();
int nnTot = 3*numNodes;
@@ -637,14 +644,13 @@ void GRbf::RbfLapSurface_global_CPM_high_2(const fullMatrix<double> &cntrs,
Msg::Info("*** RBF builded Laplacian operator");
}
-
//NEW METHOD #2 CPM GLOBAL HIGH
//Produces a nxn differentiation matrix (like the projection method)
//So the local method for this is the local projection method
void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
- fullMatrix<double> &Oper){
-
+ fullMatrix<double> &Oper)
+{
Msg::Info("*** RBF ... building Laplacian operator");
int numNodes = cntrs.size1();
int nnTot = 3*numNodes;
@@ -702,8 +708,8 @@ void GRbf::RbfLapSurface_global_CPM_high(const fullMatrix<double> &cntrs,
void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
- fullMatrix<double> &Oper) {
-
+ fullMatrix<double> &Oper)
+{
int numNodes = cntrs.size1();
int nnTot = 3*numNodes;
Oper.resize(nnTot,nnTot);
@@ -758,7 +764,8 @@ void GRbf::RbfLapSurface_global_CPM_low(const fullMatrix<double> &cntrs,
void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
std::vector<MVertex*> bcNodes,
std::vector<double> coords,
- std::map<MVertex*, SPoint3> &rbf_param){
+ std::map<MVertex*, SPoint3> &rbf_param)
+{
Msg::Info("*** RBF ... solving map");
int nb = Oper.size2();
UV.resize(nb,2);
@@ -852,13 +859,12 @@ void GRbf::solveHarmonicMap(fullMatrix<double> Oper,
Msg::Info("*** RBF solved map");
exportParametrizedMesh(UV, nbNodes);
-
}
bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
double &XX, double &YY, double &ZZ,
- SVector3 &dXdu, SVector3& dXdv, int num_neighbours){
-
+ SVector3 &dXdu, SVector3& dXdv, int num_neighbours)
+{
if (u_eval == lastU && v_eval == lastV){
XX = lastX;
YY = lastY;
@@ -875,16 +881,16 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
u_vec_eval(0,1) = v_eval;
u_vec_eval(0,2) = 0.0;
- double dist_min = 1.e6;
-
+ double dist_min = 1.e6;
+
#if defined (HAVE_ANN)
- double uvw[3] = { u_eval, v_eval, 0.0 };
- ANNidxArray index = new ANNidx[num_neighbours];
- ANNdistArray dist = new ANNdist[num_neighbours];
- UVkdtree->annkSearch(uvw, num_neighbours, index, dist);
-
- for (int i = 0; i < num_neighbours; i++){
-
+ double uvw[3] = { u_eval, v_eval, 0.0 };
+ ANNidxArray index = new ANNidx[num_neighbours];
+ ANNdistArray dist = new ANNdist[num_neighbours];
+ UVkdtree->annkSearch(uvw, num_neighbours, index, dist);
+
+ for (int i = 0; i < num_neighbours; i++){
+
u_vec(i,0) = UV(index[i],0);
u_vec(i,1) = UV(index[i],1);
u_vec(i,2) = 0.0;
@@ -898,34 +904,33 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
(UV(index[i],1)-UV(index[j],1))*(UV(index[i],1)-UV(index[j],1)));
if (dist<dist_min && dist > 1.e-5) dist_min = dist;
}
- }
- delete [] index;
- delete [] dist;
+ }
+ delete [] index;
+ delete [] dist;
#endif
- _inUV = 1;
- deltaUV = 0.3*dist_min;
- evalRbfDer(0, u_vec, u_vec_eval,xyz_local, nodes_eval);
- evalRbfDer(1, u_vec, u_vec_eval,xyz_local, xu);
- evalRbfDer(2, u_vec, u_vec_eval,xyz_local, xv);
- _inUV= 0;
-
- XX = nodes_eval(0,0)*sBox;
- YY = nodes_eval(0,1)*sBox;
- ZZ = nodes_eval(0,2)*sBox;
- dXdu = SVector3(xu(0,0)*sBox, xu(0,1)*sBox, xu(0,2)*sBox);
- dXdv = SVector3(xv(0,0)*sBox, xv(0,1)*sBox, xv(0,2)*sBox);
-
- //store last computation
- lastU = u_eval;
- lastV = v_eval;
- lastX = XX;
- lastY = YY;
- lastZ = ZZ;
- lastDXDU = dXdu ;
- lastDXDV = dXdv ;
- return true;
-
+ _inUV = 1;
+ deltaUV = 0.3*dist_min;
+ evalRbfDer(0, u_vec, u_vec_eval,xyz_local, nodes_eval);
+ evalRbfDer(1, u_vec, u_vec_eval,xyz_local, xu);
+ evalRbfDer(2, u_vec, u_vec_eval,xyz_local, xv);
+ _inUV= 0;
+
+ XX = nodes_eval(0,0)*sBox;
+ YY = nodes_eval(0,1)*sBox;
+ ZZ = nodes_eval(0,2)*sBox;
+ dXdu = SVector3(xu(0,0)*sBox, xu(0,1)*sBox, xu(0,2)*sBox);
+ dXdv = SVector3(xv(0,0)*sBox, xv(0,1)*sBox, xv(0,2)*sBox);
+
+ //store last computation
+ lastU = u_eval;
+ lastV = v_eval;
+ lastX = XX;
+ lastY = YY;
+ lastZ = ZZ;
+ lastDXDU = dXdu ;
+ lastDXDV = dXdv ;
+ return true;
}
// bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
@@ -1083,4 +1088,3 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
// #endif
// }
-
diff --git a/Geo/GRbf.h b/Geo/GRbf.h
index 6d2f89c23acabe6f606b72d9eb2cbfd7d9d567a9..bd9a01a99bb4cd2a96da28e81d400575b6d8ec00 100644
--- a/Geo/GRbf.h
+++ b/Geo/GRbf.h
@@ -1,5 +1,12 @@
-#ifndef _RBF_H
-#define _RBF_H
+// Gmsh - Copyright (C) 1997-2011 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to <gmsh@geuz.org>.
+//
+// Contributed by FIXME
+
+#ifndef _RBF_H_
+#define _RBF_H_
#include <math.h>
#include <vector>
@@ -8,11 +15,11 @@
#include "SVector3.h"
#include "MVertex.h"
#include "Context.h"
+
#if defined(HAVE_ANN)
class ANNkd_tree;
#endif
-
class Sphere{
public:
int index;
@@ -21,7 +28,6 @@ class Sphere{
};
class GRbf {
-
std::map<MVertex*, int> _mapV;
std::map<MVertex*, int> _mapAllV;
std::map<int, std::vector<int> > nodesInSphere;
@@ -70,8 +76,10 @@ class GRbf {
void buildOctree(double radius);
void buildXYZkdtree();
- // Sets up the surface generation problem as suggested by Beatson et al. Introduction of 2 extra points per node. The
- // function values at the nodes on the surface are set to 0 while the function values inside are set to -1 and outside to 1.
+ // Sets up the surface generation problem as suggested by Beatson et
+ // al. Introduction of 2 extra points per node. The function values
+ // at the nodes on the surface are set to 0 while the function
+ // values inside are set to -1 and outside to 1.
void setup_level_set(const fullMatrix<double> &cntrs,
const fullMatrix<double> &normals,
fullMatrix<double> &level_set_nodes,
@@ -81,20 +89,22 @@ class GRbf {
// This is to avoid the introduction of removable singularities at r=0.
double evalRadialFnDer (int p, double dx, double dy, double dz, double ep);
-
// Generates the RBF collocation matrix for data in d-dimensions, associated with the
//(p)th derivative of the radial function w.r.t. the (q)th variable
fullMatrix<double> generateRbfMat(int p,
const fullMatrix<double> &nodes1,
const fullMatrix<double> &nodes2);
- // Computes the interpolation(p==0) or the derivative (p!=0) operator(mxn) (n:number of centers, m: number of evaluation nodes)
+ // Computes the interpolation(p==0) or the derivative (p!=0)
+ // operator(mxn) (n:number of centers, m: number of evaluation
+ // nodes)
void RbfOp(int p, // (p)th derivatives
const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes,
fullMatrix<double> &D);
- // Computes the interpolant(p==0) or the derivative (p!=0) of the function values entered and evaluates it at the new nodes
+ // Computes the interpolant(p==0) or the derivative (p!=0) of the
+ // function values entered and evaluates it at the new nodes
void evalRbfDer(int p, // (p)th derivatives
const fullMatrix<double> &cntrs,
const fullMatrix<double> &nodes,
@@ -133,25 +143,23 @@ class GRbf {
fullMatrix<double> &D);
// Calculates the curvature of a surface at centers
-
void curvatureRBF(const fullMatrix<double> &cntrs,
fullMatrix<double> &curvature);
-
void computeCurvature(const fullMatrix<double> &cntrs,
std::map<MVertex*, double>&rbf_curv);
-
void computeLocalCurvature(const fullMatrix<double> &cntrs,
std::map<MVertex*, double>&rbf_curv);
- //Finds the U,V,S (in the 0-level set) that are the 'num_neighbours' closest to u_eval and v_eval.
- //Thus in total, we're working with '3*num_neighbours' nodes
- //Say that the vector 'index' gives us the indices of the closest points
- bool UVStoXYZ(const double u_eval, const double v_eval,
- double &XX, double &YY, double &ZZ,
- SVector3 &dXdu, SVector3& dxdv, int num_neighbours=15);
-
- void solveHarmonicMap(fullMatrix<double> Oper, std::vector<MVertex*> ordered,
- std::vector<double> coords, std::map<MVertex*, SPoint3> &rbf_param);
+ //Finds the U,V,S (in the 0-level set) that are the 'num_neighbours'
+ //closest to u_eval and v_eval. Thus in total, we're working with
+ //'3*num_neighbours' nodes Say that the vector 'index' gives us the
+ //indices of the closest points
+ bool UVStoXYZ(const double u_eval, const double v_eval,
+ double &XX, double &YY, double &ZZ,
+ SVector3 &dXdu, SVector3& dxdv, int num_neighbours=15);
+
+ void solveHarmonicMap(fullMatrix<double> Oper, std::vector<MVertex*> ordered,
+ std::vector<double> coords, std::map<MVertex*, SPoint3> &rbf_param);
inline const fullMatrix<double> getUV() {return UV;};
inline const fullMatrix<double> getXYZ() {return centers;};
diff --git a/Geo/MVertexBoundaryLayerData.cpp b/Geo/MVertexBoundaryLayerData.cpp
index 0257cfafa4e30785beb15f9c5c764de3b2ea6354..a98bd1279b6646d7a217d7308b427e8c5d54ce0f 100644
--- a/Geo/MVertexBoundaryLayerData.cpp
+++ b/Geo/MVertexBoundaryLayerData.cpp
@@ -5,7 +5,8 @@
#include "MVertexBoundaryLayerData.h"
-std::vector<MVertex*>* MVertexBoundaryLayerData::getChildren(int i){
+std::vector<MVertex*>* MVertexBoundaryLayerData::getChildren(int i)
+{
if (i < this->children.size() && i >= 0) {
return &(children[i]);
} else {
@@ -13,7 +14,8 @@ std::vector<MVertex*>* MVertexBoundaryLayerData::getChildren(int i){
}
}
-int MVertexBoundaryLayerData::getNumChildren(int i) {
+int MVertexBoundaryLayerData::getNumChildren(int i)
+{
if (i < this->children.size() && i >= 0) {
return this->children[i].size();
} else {
@@ -21,10 +23,12 @@ int MVertexBoundaryLayerData::getNumChildren(int i) {
}
}
-int MVertexBoundaryLayerData::getNumChildrenFamilies() {
+int MVertexBoundaryLayerData::getNumChildrenFamilies()
+{
return this->children.size();
}
-void MVertexBoundaryLayerData::addChildrenFamily(std::vector<MVertex*> family) {
+void MVertexBoundaryLayerData::addChildrenFamily(std::vector<MVertex*> family)
+{
this->children.push_back(family);
}
diff --git a/Geo/MVertexBoundaryLayerData.h b/Geo/MVertexBoundaryLayerData.h
index 91ff8104fcbc1583ebad21e63ba64ebcf1262309..d4bc1d2f39cf5f12fc615e8845a993d1a636a05c 100644
--- a/Geo/MVertexBoundaryLayerData.h
+++ b/Geo/MVertexBoundaryLayerData.h
@@ -10,14 +10,12 @@
class MVertex;
-/* A simple data structure to keep track of the "children" of
- * vertices in a boundary layer meshes.
- *
- * It should be filled for each MVertex on the boundary of the
- * geometry with the vertices along the normal direction, in order.
- */
-class MVertexBoundaryLayerData {
+// A simple data structure to keep track of the "children" of vertices
+// in a boundary layer mesh. It should be filled for each MVertex on
+// the boundary of the geometry with the vertices along the normal
+// direction, in order.
+class MVertexBoundaryLayerData {
private:
std::vector<std::vector<MVertex*> > children;
@@ -32,5 +30,4 @@ class MVertexBoundaryLayerData {
void addChildrenFamily(std::vector<MVertex*> family);
};
-
#endif
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index e45571e93107f917381f5eec907c7adb75fac229..ddb3c063e77e9c46b9cc2bed95eb565908079686 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -41,7 +41,6 @@ GEdge *getOCCEdgeByNativePtr(GModel *model, TopoDS_Edge toFind)
return 0;
}
-
OCCEdge::OCCEdge(GModel *model, TopoDS_Edge edge, int num, GVertex *v1, GVertex *v2)
: GEdge(model, num, v1, v2), c(edge), trimmed(0)
{
@@ -295,8 +294,8 @@ void OCCEdge::writeGEO(FILE *fp)
// in gluing operations for example
void OCCEdge::replaceEndingPointsInternals(GVertex *g0, GVertex *g1)
{
- TopoDS_Vertex aV1 = *((TopoDS_Vertex*)v0->getNativePtr());
- TopoDS_Vertex aV2 = *((TopoDS_Vertex*)v1->getNativePtr());
+ TopoDS_Vertex aV1 = *((TopoDS_Vertex*)v0->getNativePtr());
+ TopoDS_Vertex aV2 = *((TopoDS_Vertex*)v1->getNativePtr());
TopoDS_Vertex aVR1 = *((TopoDS_Vertex*)g0->getNativePtr());
TopoDS_Vertex aVR2 = *((TopoDS_Vertex*)g1->getNativePtr());
@@ -304,37 +303,32 @@ void OCCEdge::replaceEndingPointsInternals(GVertex *g0, GVertex *g1)
// v0,v1,g0,g1,v0->tag(),v1->tag(),g0->tag(),g1->tag(),tag());
Standard_Boolean bIsDE = BRep_Tool::Degenerated(c);
- //
+
TopoDS_Edge aEx = c;
aEx.Orientation(TopAbs_FORWARD);
Standard_Real t1=s0;
Standard_Real t2=s1;
- //
+
aVR1.Orientation(TopAbs_FORWARD);
aVR2.Orientation(TopAbs_REVERSED);
- //
+
if (bIsDE) {
Standard_Real aTol;
BRep_Builder aBB;
TopoDS_Edge E;
TopAbs_Orientation anOrE;
- //
- anOrE=c.Orientation();
- aTol=BRep_Tool::Tolerance(c);
- //
- E=aEx;
+ anOrE = c.Orientation();
+ aTol = BRep_Tool::Tolerance(c);
+ E = aEx;
E.EmptyCopy();
- //
- aBB.Add (E, aVR1);
- aBB.Add (E, aVR2);
+ aBB.Add(E, aVR1);
+ aBB.Add(E, aVR2);
aBB.Range(E, t1, t2);
aBB.Degenerated(E, Standard_True);
aBB.UpdateEdge(E, aTol);
- //
_replacement=E;
}
- //
else {
BOPTools_Tools::MakeSplitEdge(aEx, aVR1, t1, aVR2, t2, _replacement);
}
@@ -345,7 +339,6 @@ void OCCEdge::replaceEndingPointsInternals(GVertex *g0, GVertex *g1)
//build the reverse curve
c_rev = c;
c_rev.Reverse();
-
}
#endif
diff --git a/Geo/OCCEdge.h b/Geo/OCCEdge.h
index a555434d41928fc81bb2a9e25618f983e3dae668..cebd24038d0e70e047e3355ae19ac98cfdb77342 100644
--- a/Geo/OCCEdge.h
+++ b/Geo/OCCEdge.h
@@ -48,6 +48,7 @@ class OCCEdge : public GEdge {
TopoDS_Edge getTopoDS_EdgeOld() const {return _replacement;}
void replaceEndingPointsInternals(GVertex *, GVertex *);
};
+
GEdge *getOCCEdgeByNativePtr(GModel *model, TopoDS_Edge toFind);
#endif
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index ebf4db718a2022a9ec4b4df38ed4323756901263..7ea79f57be2c6389ba61125deca28df497e90035 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -53,7 +53,9 @@ void OCCFace::setup()
}
else{
l_wire.push_back(e);
- Msg::Debug("Edge %d (%d --> %d) ori %d", e->tag(), e->getBeginVertex()->tag(), e->getEndVertex()->tag(), edge.Orientation());
+ Msg::Debug("Edge %d (%d --> %d) ori %d", e->tag(),
+ e->getBeginVertex()->tag(), e->getEndVertex()->tag(),
+ edge.Orientation());
e->addFace(this);
if(!e->is3D()){
OCCEdge *occe = (OCCEdge*)e;
@@ -63,7 +65,7 @@ void OCCFace::setup()
}
GEdgeLoop el(l_wire);
- // printf("l_wire of size %d %d\n",l_wire.size(),el.count());
+ // printf("l_wire of size %d %d\n",l_wire.size(),el.count());
for(GEdgeLoop::citer it = el.begin(); it != el.end(); ++it){
l_edges.push_back(it->ge);
l_dirs.push_back(it->_sign);
@@ -96,11 +98,12 @@ void OCCFace::setup()
vmax += fabs(dv) / 100.0;
occface = BRep_Tool::Surface(s);
- // std::list<GEdge*>::const_iterator it = l_edges.begin();
- // for (; it != l_edges.end(); ++it){
- // printf("edge %d : %d %d iseam %d \n",(*it)->tag(),(*it)->getBeginVertex()->tag(),(*it)->getEndVertex()->tag(),(*it)->isSeam(this));
- // }
-
+ // std::list<GEdge*>::const_iterator it = l_edges.begin();
+ // for (; it != l_edges.end(); ++it){
+ // printf("edge %d : %d %d iseam %d \n", (*it)->tag(),
+ // (*it)->getBeginVertex()->tag(), (*it)->getEndVertex()->tag(),
+ // (*it)->isSeam(this));
+ // }
}
Range<double> OCCFace::parBounds(int i) const
@@ -133,7 +136,6 @@ Pair<SVector3,SVector3> OCCFace::firstDer(const SPoint2 ¶m) const
gp_Pnt pnt;
gp_Vec du, dv;
occface->D1(param.x(), param.y(), pnt, du, dv);
-
return Pair<SVector3,SVector3>(SVector3(du.X(), du.Y(), du.Z()),
SVector3(dv.X(), dv.Y(), dv.Z()));
}
@@ -170,7 +172,7 @@ GPoint OCCFace::closestPoint(const SPoint3 &qp, const double initialGuess[2]) co
double pp[2] = {initialGuess[0],initialGuess[1]};
proj.LowerDistanceParameters(pp[0], pp[1]);
- // Msg::Info("projection lower distance parameters %g %g",pp[0],pp[1]);
+ // Msg::Info("projection lower distance parameters %g %g",pp[0],pp[1]);
if((pp[0] < umin || umax < pp[0]) || (pp[1]<vmin || vmax<pp[1])){
Msg::Error("Point projection is out of face bounds");
@@ -217,18 +219,16 @@ GEntity::GeomType OCCFace::geomType() const
double OCCFace::curvatureMax(const SPoint2 ¶m) const
{
-
const double eps = 1.e-12;
BRepAdaptor_Surface sf(s, Standard_True);
BRepLProp_SLProps prop(sf, 2, eps);
- prop.SetParameters (param.x(),param.y());
+ prop.SetParameters (param.x(), param.y());
- if (!prop.IsCurvatureDefined()){
+ if(!prop.IsCurvatureDefined()){
return eps;
}
return std::max(fabs(prop.MinCurvature()), fabs(prop.MaxCurvature()));
-
}
double OCCFace::curvatures(const SPoint2 ¶m,
@@ -506,7 +506,6 @@ void OCCFace::replaceEdgesInternal(std::list<GEdge*> &new_edges)
else {
aER.Orientation(aE.Orientation());
}
- //
aBB.Add(newWire, aER);
}
aBB.Add(newFace, newWire);
diff --git a/Geo/OCCFace.h b/Geo/OCCFace.h
index 9fa9ffbe8f10402d74c9c2f1a1cb5f8235cb7c6b..752735c27ba4ea93d7a31abeeaf970c7ab8ac05c 100644
--- a/Geo/OCCFace.h
+++ b/Geo/OCCFace.h
@@ -52,6 +52,7 @@ class OCCFace : public GFace {
// tells if it's a sphere, and if it is, returns parameters
virtual bool isSphere (double &radius, SPoint3 ¢er) const;
};
+
GFace *getOCCFaceByNativePtr(GModel *model, TopoDS_Face toFind);
#endif
diff --git a/Geo/OCCRegion.h b/Geo/OCCRegion.h
index 8be01f5eef681621f13b3426fce5db899c7f1413..211cb301b79b722041f50f58a2309b883df71f7b 100644
--- a/Geo/OCCRegion.h
+++ b/Geo/OCCRegion.h
@@ -26,6 +26,7 @@ class OCCRegion : public GRegion {
};
GRegion *getOCCRegionByNativePtr(GModel *model, TopoDS_Solid toFind);
+
#endif
#endif
diff --git a/Geo/OCCVertex.h b/Geo/OCCVertex.h
index 1ff76ef359515e669507e89a4967b7614757b488..b0d4d4b19f87029a538e610f2f2a8e4153052ecb 100644
--- a/Geo/OCCVertex.h
+++ b/Geo/OCCVertex.h
@@ -32,6 +32,7 @@ class OCCVertex : public GVertex {
virtual SPoint2 reparamOnFace(const GFace *gf, int) const;
TopoDS_Vertex getShape() { return v; }
};
+
GVertex *getOCCVertexByNativePtr(GModel *model, TopoDS_Vertex toFind);
#endif
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index de5974559e7fe7415677e67cb7c93bd0edce91d7..cb3d339157033dfb297a9d23e7eefeca440933dc 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -15,6 +15,7 @@
#include <stdlib.h> // for abs()
#include <vector>
#include "fullMatrix.h"
+#include "GModel.h"
#include "MVertex.h"
#include "GmshConfig.h"
#include "mathEvaluator.h"
@@ -254,7 +255,6 @@ public:
int type() const { return UNKNOWN; }
};
-class GModel;
class gLevelsetDistGeom: public gLevelsetPrimitive
{
GModel *_model;
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 19fc0c1c0b09656457e222073119db577397dd79..5d142c8f35b872007e46039826bcc4b9f0b0adc9 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -15,7 +15,7 @@
void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace *gf)
{
- FILE *f = fopen(iii, "w");
+ FILE *f = fopen(iii, "w");
if(!f) return;
fprintf(f, "View \"scalar\" {\n");
std::list<BDS_Face*>::iterator tit = t.begin();
@@ -263,7 +263,8 @@ BDS_Edge *BDS_Mesh::recover_edge_fast(BDS_Point *p1, BDS_Point *p2){
return 0;
}
-BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal, std::set<EdgeToRecover> *e2r,
+BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal,
+ std::set<EdgeToRecover> *e2r,
std::set<EdgeToRecover> *not_recovered)
{
BDS_Edge *e = find_edge(num1, num2);
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index c2c84d6180f14b9bf95d6043459443a37edd0583..79990250de8175f12e6b9a18553b89391b1076a6 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -30,7 +30,6 @@ static void computeMeshMetricsForBamg(GFace *gf, int numV,
double *mm11, double *mm12, double *mm22,
int iter)
{
-
// char name[245];
// sprintf(name,"bgmBamg-%d-%d.pos",gf->tag(),iter);
// if (iter < 2){
@@ -67,7 +66,6 @@ static void computeMeshMetricsForBamg(GFace *gf, int numV,
static void meshGFaceBamg_(GFace *gf, int iter, bool initialMesh)
{
-
std::set<MVertex*> all;
std::map<int,MVertex*> recover;
for (unsigned int i = 0; i < gf->triangles.size(); i++){
@@ -215,8 +213,8 @@ void meshGFaceBamg(GFace *gf)
void meshGFaceBamg(GFace *gf)
{
- Msg::Error("Gmsh msust be compiled with the Bidimensional Anisotropic Mesh "
- "Generator (Bamg) in order to support that option");
+ Msg::Error("This version of Gmsh is not compiled with Bidimensional "
+ "Anisotropic Mesh Generator (BAMG) support");
}
#endif
diff --git a/Mesh/meshGFaceBamg.h b/Mesh/meshGFaceBamg.h
index 2bff95afe227bfae357ad22a4f712dbe558186e4..d786b1aa5fdc22394a987f5f5e131e38b9e4c6b8 100644
--- a/Mesh/meshGFaceBamg.h
+++ b/Mesh/meshGFaceBamg.h
@@ -7,6 +7,7 @@
#define _MESHGFACE_BAMG_
class GFace;
+
void meshGFaceBamg(GFace *gf);
#endif
diff --git a/Mesh/meshGRegionMMG3D.cpp b/Mesh/meshGRegionMMG3D.cpp
index 6cda27bc004808876b7eae735f0ec2f3dde5e58c..fa1f440a24250f25afe2252eb4d762a0ace984fc 100644
--- a/Mesh/meshGRegionMMG3D.cpp
+++ b/Mesh/meshGRegionMMG3D.cpp
@@ -23,7 +23,7 @@ extern "C" {
#define M_UNUSED (1 << 0)
}
-void MMG2gmsh (GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gmsh)
+static void MMG2gmsh(GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gmsh)
{
std::map<int,MVertex*> kToMVertex;
for (int k=1;k<= mmg->np ; k++){
@@ -37,7 +37,6 @@ void MMG2gmsh (GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gmsh)
}
else kToMVertex[k] = it->second;
}
-
for (int k=1; k<=mmg->ne; k++) {
MMG_pTetra ptetra = &mmg->tetra[k];
@@ -55,7 +54,8 @@ void MMG2gmsh (GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gmsh)
}
}
-void gmsh2MMG (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*> &mmg2gmsh)
+static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
+ std::map<int,MVertex*> &mmg2gmsh)
{
mmg->ne = gr->tetrahedra.size();
std::set<MVertex*> allVertices;
@@ -112,7 +112,8 @@ void gmsh2MMG (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*>
int k=1;
int count = 1;//sol->offset;
std::map<int,int> gmsh2mmg_num;
- for (std::set<MVertex*>::iterator it = allVertices.begin() ; it != allVertices.end() ; ++it){
+ for (std::set<MVertex*>::iterator it = allVertices.begin();
+ it != allVertices.end(); ++it){
MMG_pPoint ppt = &mmg->point[k];
ppt->c[0] = (*it)->x();
@@ -184,7 +185,7 @@ void gmsh2MMG (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*>
}
-void updateSizes (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol)
+static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol)
{
std::list<GFace*> f = gr->faces();
/*
@@ -238,8 +239,7 @@ void updateSizes (GRegion *gr, MMG_pMesh mmg, MMG_pSol sol)
sol->metold = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
}
-
-void FREEMMG (MMG_pMesh mmgMesh, MMG_pSol mmgSol)
+static void freeMMG(MMG_pMesh mmgMesh, MMG_pSol mmgSol)
{
free(mmgMesh->point);
// free(mmgMesh->disp->alpha);
@@ -275,18 +275,18 @@ void refineMeshMMG(GRegion *gr)
gr->deleteVertexArrays();
for (int i=0;i<gr->tetrahedra.size();++i)delete gr->tetrahedra[i];
gr->tetrahedra.clear();
- // for (int i=0;i<gr->mesh_vertices.size();++i)delete gr->mesh_vertices[i];
- //gr->mesh_vertices.clear();
+ // for (int i=0;i<gr->mesh_vertices.size();++i)delete gr->mesh_vertices[i];
+ // gr->mesh_vertices.clear();
- MMG2gmsh (gr, mmg, mmg2gmsh);
- FREEMMG (mmg,sol);
+ MMG2gmsh(gr, mmg, mmg2gmsh);
+ freeMMG(mmg, sol);
}
#else
void refineMeshMMG(GRegion *gr)
{
- Msg::Error("You should compile your version of Gmsh with MMG3D, the Mobile Mesh Generator");
+ Msg::Error("This version of Gmsh is not compiled with MMG3D support");
}
#endif
diff --git a/Mesh/meshGRegionMMG3D.h b/Mesh/meshGRegionMMG3D.h
index 2567534d2444a48d780161649047acacb30cde17..698ad5053f323ff85ecb1f1811e333b023b61f24 100644
--- a/Mesh/meshGRegionMMG3D.h
+++ b/Mesh/meshGRegionMMG3D.h
@@ -5,6 +5,9 @@
#ifndef _MESHGREGIONMMG3D_H_
#define _MESHGREGIONMMG3D_H_
+
class GRegion;
+
void refineMeshMMG(GRegion *gr);
+
#endif
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index 4204803abf262912820aef56260c58956bf7daa7..7ba924e10aff64906035f5f71b9c3574640bd44f 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -21,10 +21,10 @@
#include "robustPredicates.h"
#include "BackgroundMesh.h"
-// FIXME: should not introduce dependencies to Geo/ code in Numeric
+// FIXME: should not introduce dependencies to Geo/ code in Numeric!
#include "GPoint.h"
#include "GFace.h"
-#include "GFaceCompound.h"
+#include "GEdgeCompound.h"
#include "MLine.h"
#define Pred(x) ((x)->prev)
@@ -551,8 +551,7 @@ void DocRecord::voronoiCell(PointNumero pt, std::vector<SPoint2> &pts) const
}
/*
- Consider N points {X_1, \dots, X_N}
- We want to find the point X_P that verifies
+ Consider N points {X_1, \dots, X_N}. We want to find the point X_P that verifies
min max | X_i - X_P | , i=1,\dots,N
@@ -567,7 +566,6 @@ void DocRecord::voronoiCell(PointNumero pt, std::vector<SPoint2> &pts) const
=> 2 \int_V || X - X_P||^{2m-1} dv = 0 => X_P is somewhere ...
now in infinite norm, how to find X_P ?
-
*/
void DocRecord::makePosView(std::string fileName, GFace *gf)
@@ -866,61 +864,69 @@ void DocRecord::setPoints(fullMatrix<double> *p)
}
}
-void DocRecord::recur_tag_triangles(int iTriangle,std::set<int>& taggedTriangles,std::map<std::pair<void*,void*>,std::pair<int,int> >& edgesToTriangles)
+void DocRecord::recur_tag_triangles(int iTriangle, std::set<int>& taggedTriangles,
+ std::map<std::pair<void*, void*>,
+ std::pair<int,int> > &edgesToTriangles)
{
- if(taggedTriangles.find(iTriangle)!=taggedTriangles.end()) return;
+ if(taggedTriangles.find(iTriangle) != taggedTriangles.end()) return;
taggedTriangles.insert(iTriangle);
-
+
int a = triangles[iTriangle].a;
int b = triangles[iTriangle].b;
int c = triangles[iTriangle].c;
- PointRecord* p[3] = {&points[a],&points[b],&points[c]};
+ PointRecord* p[3] = {&points[a], &points[b], &points[c]};
for(int j=0;j<3;j++){
- if(!lookForBoundaryEdge(p[j]->data,p[(j+1)%3]->data)){
- std::pair<void*,void*> ab = std::make_pair(std::min(p[j]->data,p[(j+1)%3]->data),std::max(p[j]->data,p[(j+1)%3]->data));
- std::pair<int,int>& adj = (edgesToTriangles.find(ab))->second;
- if(iTriangle==adj.first && adj.second!=-1) recur_tag_triangles(adj.second,taggedTriangles,edgesToTriangles);
- else recur_tag_triangles(adj.first,taggedTriangles,edgesToTriangles);
- }
+ if(!lookForBoundaryEdge(p[j]->data, p[(j + 1) % 3]->data)){
+ std::pair<void*,void*> ab = std::make_pair
+ (std::min(p[j]->data, p[(j + 1) % 3]->data),
+ std::max(p[j]->data, p[(j + 1) % 3]->data));
+ std::pair<int,int>& adj = (edgesToTriangles.find(ab))->second;
+ if(iTriangle == adj.first && adj.second != -1)
+ recur_tag_triangles(adj.second, taggedTriangles, edgesToTriangles);
+ else
+ recur_tag_triangles(adj.first, taggedTriangles, edgesToTriangles);
+ }
}
}
-std::set<int> DocRecord::tagInterior(double x,double y)
+std::set<int> DocRecord::tagInterior(double x, double y)
{
std::map<std::pair<void*,void*>,std::pair<int,int> > edgesToTriangles;
int iFirst = 1;
- for(int i=0;i<numTriangles;i++){
+ for(int i = 0; i < numTriangles; i++){
int a = triangles[i].a;
- int b = triangles[i].b;
- int c = triangles[i].c;
- PointRecord* p[3] = {&points[a],&points[b],&points[c]};
-
- //Weisstein, Eric W. "Triangle Interior." From MathWorld--A Wolfram Web Resource.
- double x0 = points[a].where.h;
- double y0 = points[a].where.v;
- double x1 = points[b].where.h - points[a].where.h;
- double y1 = points[b].where.v - points[a].where.v;
- double x2 = points[c].where.h - points[a].where.h;
- double y2 = points[c].where.v - points[a].where.v;
- double k1 = ((x*y2-x2*y)-(x0*y2-x2*y0))/(x1*y2-x2*y1);
- double k2 = -((x*y1-x1*y)-(x0*y1-x1*y0))/(x1*y2-x2*y1);
- if(k1>0.0 && k2>0.0 && k1+k2<1.0) iFirst = i;
-
- for(int j=0;j<3;j++){
- std::pair<void*,void*> ab = std::make_pair(std::min(p[j]->data,p[(j+1)%3]->data),std::max(p[j]->data,p[(j+1)%3]->data));
- std::map<std::pair<void*,void*>,std::pair<int,int> >::iterator it = edgesToTriangles.find(ab);
- if(it==edgesToTriangles.end()){
- edgesToTriangles[ab] = std::make_pair(i,-1);
- }
- else{
- it->second.second = i;
- }
- }
+ int b = triangles[i].b;
+ int c = triangles[i].c;
+ PointRecord* p[3] = {&points[a],&points[b],&points[c]};
+
+ //Weisstein, Eric W. "Triangle Interior." From MathWorld--A Wolfram Web Resource.
+ double x0 = points[a].where.h;
+ double y0 = points[a].where.v;
+ double x1 = points[b].where.h - points[a].where.h;
+ double y1 = points[b].where.v - points[a].where.v;
+ double x2 = points[c].where.h - points[a].where.h;
+ double y2 = points[c].where.v - points[a].where.v;
+ double k1 = ((x*y2-x2*y)-(x0*y2-x2*y0))/(x1*y2-x2*y1);
+ double k2 = -((x*y1-x1*y)-(x0*y1-x1*y0))/(x1*y2-x2*y1);
+ if(k1>0.0 && k2>0.0 && k1+k2<1.0) iFirst = i;
+ for(int j=0;j<3;j++){
+ std::pair<void*,void*> ab = std::make_pair
+ (std::min(p[j]->data, p[(j + 1) % 3]->data),
+ std::max(p[j]->data, p[(j + 1) % 3]->data));
+ std::map<std::pair<void*,void*>, std::pair<int,int> >::iterator it =
+ edgesToTriangles.find(ab);
+ if(it == edgesToTriangles.end()){
+ edgesToTriangles[ab] = std::make_pair(i, -1);
+ }
+ else{
+ it->second.second = i;
+ }
+ }
}
std::set<int> taggedTriangles;
- recur_tag_triangles(iFirst,taggedTriangles,edgesToTriangles);
+ recur_tag_triangles(iFirst, taggedTriangles, edgesToTriangles);
return taggedTriangles;
}
@@ -941,35 +947,32 @@ void DocRecord::concave(double x,double y,GFace* gf)
list = gf->edges();
replaceMeshCompound(gf,list);
- for(it1=list.begin();it1!=list.end();it1++){
- edge = *it1;
- for(i=0;i<edge->getNumMeshElements();i++){
- element = edge->getMeshElement(i);
- vertex1 = element->getVertex(0);
- vertex2 = element->getVertex(1);
- addBoundaryEdge(vertex1,vertex2);
- }
+ for(it1 = list.begin(); it1 != list.end(); it1++){
+ edge = *it1;
+ for(i=0;i<edge->getNumMeshElements();i++){
+ element = edge->getMeshElement(i);
+ vertex1 = element->getVertex(0);
+ vertex2 = element->getVertex(1);
+ addBoundaryEdge(vertex1,vertex2);
+ }
}
- for(i=0;i<numPoints;i++){
+ for(i = 0;i < numPoints; i++){
points[i].vicinity.clear();
}
MakeMeshWithPoints();
set = tagInterior(x,y);
- for(it2=set.begin();it2!=set.end();it2++){
- index1 = triangles[*it2].a;
- index2 = triangles[*it2].b;
- index3 = triangles[*it2].c;
-
- add(index1,index2);
- add(index1,index3);
-
- add(index2,index1);
- add(index2,index3);
-
- add(index3,index1);
- add(index3,index2);
+ for(it2 = set.begin(); it2 != set.end(); it2++){
+ index1 = triangles[*it2].a;
+ index2 = triangles[*it2].b;
+ index3 = triangles[*it2].c;
+ add(index1,index2);
+ add(index1,index3);
+ add(index2,index1);
+ add(index2,index3);
+ add(index3,index1);
+ add(index3,index2);
}
}
@@ -980,13 +983,15 @@ bool DocRecord::contain(int index1,int index2,int index3)
void* dataB;
dataA = points[index2].data;
dataB = points[index3].data;
- for(i=0;i<points[index1].vicinity.size()-1;i=i+2){
- if(points[index1].vicinity[i]==dataA && points[index1].vicinity[i+1]==dataB){
- return 1;
- }
- else if(points[index1].vicinity[i]==dataB && points[index1].vicinity[i+1]==dataA){
- return 1;
- }
+ for(i = 0; i < points[index1].vicinity.size() - 1; i = i+2){
+ if(points[index1].vicinity[i] == dataA &&
+ points[index1].vicinity[i + 1] == dataB){
+ return 1;
+ }
+ else if(points[index1].vicinity[i] == dataB &&
+ points[index1].vicinity[i + 1] == dataA){
+ return 1;
+ }
}
return 0;
}
@@ -1001,8 +1006,8 @@ void DocRecord::add(int index1,int index2)
void DocRecord::initialize()
{
int i;
- for(i=0;i<numPoints;i++){
- points[i].flag = 0;
+ for(i = 0; i < numPoints; i++){
+ points[i].flag = 0;
}
}
@@ -1022,22 +1027,22 @@ void DocRecord::remove_all()
int numPoints2;
PointRecord* points2;
numPoints2 = 0;
- for(i=0;i<numPoints;i++){
- if(points[i].flag==0){
+ for(i = 0; i < numPoints; i++){
+ if(points[i].flag == 0){
numPoints2++;
}
}
points2 = new PointRecord[numPoints2];
index = 0;
- for(i=0;i<numPoints;i++){
- if(points[i].flag==0){
- points2[index].where.h = points[i].where.h;
- points2[index].where.v = points[i].where.v;
- points2[index].data = points[i].data;
- points2[index].flag = points[i].flag;
- points2[index].identificator = points[i].identificator;
- index++;
- }
+ for(i = 0; i < numPoints; i++){
+ if(points[i].flag == 0){
+ points2[index].where.h = points[i].where.h;
+ points2[index].where.v = points[i].where.v;
+ points2[index].data = points[i].data;
+ points2[index].flag = points[i].flag;
+ points2[index].identificator = points[i].identificator;
+ index++;
+ }
}
delete [] points;
points = points2;
@@ -1049,9 +1054,9 @@ void DocRecord::add_point(double x,double y,GFace*face)
PointRecord point;
point.where.h = x;
point.where.v = y;
- point.data = new MVertex(x,y,0.0,(GEntity*)face,2);
+ point.data = new MVertex(x, y, 0.0, (GEntity*)face, 2);
points[numPoints] = point;
- numPoints = numPoints+1;
+ numPoints = numPoints + 1;
}
void DocRecord::build_edges()
@@ -1063,14 +1068,14 @@ void DocRecord::build_edges()
MVertex* vertex1;
MVertex* vertex2;
- for(i=0;i<numPoints;i++){
+ for( i = 0; i < numPoints; i++){
vertex1 = (MVertex*)points[i].data;
- num = _adjacencies[i].t_length;
- for(j=0;j<num;j++){
- index = _adjacencies[i].t[j];
- vertex2 = (MVertex*)points[index].data;
- add_edge(vertex1,vertex2);
- }
+ num = _adjacencies[i].t_length;
+ for(j = 0; j < num; j++){
+ index = _adjacencies[i].t[j];
+ vertex2 = (MVertex*)points[index].data;
+ add_edge(vertex1, vertex2);
+ }
}
}
@@ -1092,15 +1097,15 @@ bool DocRecord::delaunay_conformity(GFace* gf)
list = gf->edges();
replaceMeshCompound(gf,list);
- for(it=list.begin();it!=list.end();it++){
+ for(it = list.begin(); it!= list.end(); it++){
edge = *it;
- for(i=0;i<edge->getNumMeshElements();i++){
- element = edge->getMeshElement(i);
- vertex1 = element->getVertex(0);
- vertex2 = element->getVertex(1);
- flag = find_edge(vertex1,vertex2);
- if(!flag) return false;
- }
+ for(i = 0; i < edge->getNumMeshElements(); i++){
+ element = edge->getMeshElement(i);
+ vertex1 = element->getVertex(0);
+ vertex2 = element->getVertex(1);
+ flag = find_edge(vertex1,vertex2);
+ if(!flag) return false;
+ }
}
return true;
}
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index c26ff04f4a9cb42e018aa8ffe8cae0468993879d..1907f2dcf18e3a96f462a7eb8c0261e362a9122b 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -10,12 +10,11 @@
#include <algorithm>
#include "fullMatrix.h"
-// FIXME: should not introduce dependencies to Geo/ code in Numeric
+// FIXME: should not introduce dependencies to Geo/ code in Numeric!
#include "SPoint2.h"
#include "simpleFunction.h"
#include "Octree.h"
-class binding;
class GFace;
typedef struct _CDLIST DListRecord, *DListPeek;
typedef int PointNumero;
@@ -108,25 +107,28 @@ class DocRecord{
std::set<std::pair<void*,void*> > boundaryEdges;
- void addBoundaryEdge(void* p1,void* p2){
- void* a = (p1 < p2) ? p1 : p2;
- void* b = (p1 > p2) ? p1 : p2;
- boundaryEdges.insert(std::make_pair(a,b));
+ void addBoundaryEdge(void* p1,void* p2)
+ {
+ void *a = (p1 < p2) ? p1 : p2;
+ void *b = (p1 > p2) ? p1 : p2;
+ boundaryEdges.insert(std::make_pair(a,b));
}
-
- bool lookForBoundaryEdge(void* p1,void* p2){
- void* a = (p1 < p2) ? p1 : p2;
- void* b = (p1 > p2) ? p1 : p2;
- std::set<std::pair<void*,void*> >::iterator it = boundaryEdges.find(std::make_pair(a,b));
- return it!=boundaryEdges.end();
+ bool lookForBoundaryEdge(void* p1,void* p2)
+ {
+ void *a = (p1 < p2) ? p1 : p2;
+ void *b = (p1 > p2) ? p1 : p2;
+ std::set<std::pair<void*,void*> >::iterator it =
+ boundaryEdges.find(std::make_pair(a,b));
+ return it != boundaryEdges.end();
}
- void recur_tag_triangles(int,std::set<int>&,std::map<std::pair<void*,void*>,std::pair<int,int> >&);
+ void recur_tag_triangles(int, std::set<int>&, std::map<std::pair<void*,void*>,
+ std::pair<int,int> >&);
std::set<int> tagInterior(double,double);
void concave(double,double,GFace*);
bool contain(int,int,int);
void add(int,int);
-
+
void initialize();
bool remove_point(int);
void remove_all();
@@ -134,17 +136,19 @@ class DocRecord{
std::set<std::pair<void*,void*> > mesh_edges;
- void add_edge(void* p1,void* p2){
- void* a = (p1 < p2) ? p1 : p2;
- void* b = (p1 > p2) ? p1 : p2;
- mesh_edges.insert(std::make_pair(a,b));
+ void add_edge(void* p1,void* p2)
+ {
+ void *a = (p1 < p2) ? p1 : p2;
+ void *b = (p1 > p2) ? p1 : p2;
+ mesh_edges.insert(std::make_pair(a,b));
}
-
- bool find_edge(void* p1,void* p2){
- void* a = (p1 < p2) ? p1 : p2;
- void* b = (p1 > p2) ? p1 : p2;
- std::set<std::pair<void*,void*> >::iterator it = mesh_edges.find(std::make_pair(a,b));
- return it!=mesh_edges.end();
+ bool find_edge(void* p1,void* p2)
+ {
+ void *a = (p1 < p2) ? p1 : p2;
+ void *b = (p1 > p2) ? p1 : p2;
+ std::set<std::pair<void*,void*> >::iterator it =
+ mesh_edges.find(std::make_pair(a,b));
+ return it != mesh_edges.end();
}
void build_edges();
@@ -168,5 +172,4 @@ void centroidOfPolygon(SPoint2 &pc,
double &areaCell,
simpleFunction<double> *bgm = 0);
-
#endif
diff --git a/Numeric/GaussIntegration.cpp b/Numeric/GaussIntegration.cpp
index 9e1f471012dc4ff8b9f20438726d4e4c19ce4caa..848ea626247973722530b66871077bfc2d413252 100644
--- a/Numeric/GaussIntegration.cpp
+++ b/Numeric/GaussIntegration.cpp
@@ -6,58 +6,65 @@ static void pts2fullMatrix(int npts, IntPt *pts, fullMatrix<double> &pMatrix,
{
pMatrix.resize(npts,3);
wMatrix.resize(npts,1);
- for (int i=0;i<npts;i++) {
- pMatrix(i,0) = pts[i].pt[0];
- pMatrix(i,1) = pts[i].pt[1];
- pMatrix(i,2) = pts[i].pt[2];
- wMatrix(i,0) = pts[i].weight;
+ for (int i = 0; i < npts; i++) {
+ pMatrix(i, 0) = pts[i].pt[0];
+ pMatrix(i, 1) = pts[i].pt[1];
+ pMatrix(i, 2) = pts[i].pt[2];
+ wMatrix(i, 0) = pts[i].weight;
}
}
-void gaussIntegration::getTriangle(int order, fullMatrix<double> &pts, fullMatrix<double> &weights)
+void gaussIntegration::getTriangle(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights)
{
pts2fullMatrix(getNGQTPts(order),getGQTPts(order),pts,weights);
}
-void gaussIntegration::getLine(int order, fullMatrix<double> &pts, fullMatrix<double> &weights)
+void gaussIntegration::getLine(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights)
{
pts2fullMatrix(getNGQLPts(order),getGQLPts(order),pts,weights);
}
-void gaussIntegration::getQuad(int order, fullMatrix<double> &pts, fullMatrix<double> &weights)
+void gaussIntegration::getQuad(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights)
{
pts2fullMatrix(getNGQQPts(order),getGQQPts(order),pts,weights);
}
-void gaussIntegration::getTetrahedron(int order, fullMatrix<double> &pts, fullMatrix<double> &weights)
+void gaussIntegration::getTetrahedron(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights)
{
pts2fullMatrix(getNGQTetPts(order),getGQTetPts(order),pts,weights);
}
-void gaussIntegration::getHexahedron(int order, fullMatrix<double> &pts, fullMatrix<double> &weights)
+void gaussIntegration::getHexahedron(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights)
{
pts2fullMatrix(getNGQHPts(order),getGQHPts(order),pts,weights);
}
-void gaussIntegration::getPrism(int order, fullMatrix<double> &pts, fullMatrix<double> &weights)
+void gaussIntegration::getPrism(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights)
{
pts2fullMatrix(getNGQPriPts(order),getGQPriPts(order),pts,weights);
}
-void gaussIntegration::get(int elementType, int order, fullMatrix<double> &pts, fullMatrix<double> &weights)
+void gaussIntegration::get(int elementType, int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights)
{
switch (elementType) {
- case TYPE_TRI : pts2fullMatrix(getNGQTPts(order), getGQTPts(order), pts, weights); break;
- case TYPE_LIN : pts2fullMatrix(getNGQLPts(order), getGQLPts(order), pts, weights); break;
- case TYPE_QUA : pts2fullMatrix(getNGQQPts(order), getGQQPts(order), pts, weights); break;
- case TYPE_TET : pts2fullMatrix(getNGQTetPts(order), getGQTetPts(order), pts, weights); break;
- case TYPE_HEX : pts2fullMatrix(getNGQHPts(order), getGQHPts(order), pts, weights); break;
- case TYPE_PRI : pts2fullMatrix(getNGQPriPts(order), getGQPriPts(order), pts, weights); break;
- case TYPE_PNT :
- weights.resize(1,1);
- weights(0,0) = 1.;
- pts.resize(1,3);
+ case TYPE_TRI : pts2fullMatrix(getNGQTPts(order), getGQTPts(order), pts, weights); break;
+ case TYPE_LIN : pts2fullMatrix(getNGQLPts(order), getGQLPts(order), pts, weights); break;
+ case TYPE_QUA : pts2fullMatrix(getNGQQPts(order), getGQQPts(order), pts, weights); break;
+ case TYPE_TET : pts2fullMatrix(getNGQTetPts(order), getGQTetPts(order), pts, weights); break;
+ case TYPE_HEX : pts2fullMatrix(getNGQHPts(order), getGQHPts(order), pts, weights); break;
+ case TYPE_PRI : pts2fullMatrix(getNGQPriPts(order), getGQPriPts(order), pts, weights); break;
+ case TYPE_PNT :
+ weights.resize(1,1);
+ weights(0,0) = 1.;
+ pts.resize(1,3);
break;
- default : Msg::Error("No integration rules defined for type %i", elementType);
+ default : Msg::Error("No integration rules defined for type %i", elementType);
}
}
diff --git a/Numeric/GaussIntegration.h b/Numeric/GaussIntegration.h
index 099b6b22af8d8734f7721f970a736f0e5c505d07..de57a0dce5f3d78dcc43078cd481b0d5ddee44f7 100644
--- a/Numeric/GaussIntegration.h
+++ b/Numeric/GaussIntegration.h
@@ -36,13 +36,20 @@ IntPt *getGQHPts(int order);
class gaussIntegration {
public:
- static void get(int elementType, int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
- static void getTriangle(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
- static void getLine(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
- static void getQuad(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
- static void getTetrahedron(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
- static void getHexahedron(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
- static void getPrism(int order, fullMatrix<double> &pts, fullMatrix<double> &weights);
+ static void get(int elementType, int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights);
+ static void getTriangle(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights);
+ static void getLine(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights);
+ static void getQuad(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights);
+ static void getTetrahedron(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights);
+ static void getHexahedron(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights);
+ static void getPrism(int order, fullMatrix<double> &pts,
+ fullMatrix<double> &weights);
};
#endif
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 6080b87e8f71286c6a6fdbde17b002324ffb7cb2..8294202c27424a380862b575ab5a21eaff0be245 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -15,9 +15,9 @@ class bezierBasis {
public :
int numLagPts;
int numDivisions;
- // The 'numLagPts' first exponents are related to 'Lagrangian' values
- fullMatrix<double> exponents; //exposants of Bezier FF
- fullMatrix<double> points; //sampling point
+ // the 'numLagPts' first exponents are related to 'Lagrangian' values
+ fullMatrix<double> exponents; // exponents of Bezier FF
+ fullMatrix<double> points; // sampling point
fullMatrix<double> matrixLag2Bez, matrixBez2Lag;
fullMatrix<double> subDivisor;
static const bezierBasis *find(int);
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index fc48c43e8c8d08ff84bbf478eb583020fd1c66ce..2051bd25287901b2e5897cfbc7ced3b0f86c8345 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -504,10 +504,10 @@ void eigenvalue(double mat[3][3], double v[3])
c[1] = 0.5 * (c[2] * c[2] - trace2(mat));
c[0] = - det3x3(mat);
- //printf("%g %g %g\n", mat[0][0], mat[0][1], mat[0][2]);
- //printf("%g %g %g\n", mat[1][0], mat[1][1], mat[1][2]);
- //printf("%g %g %g\n", mat[2][0], mat[2][1], mat[2][2]);
- //printf("%g x^3 + %g x^2 + %g x + %g = 0\n", c[3], c[2], c[1], c[0]);
+ // printf("%g %g %g\n", mat[0][0], mat[0][1], mat[0][2]);
+ // printf("%g %g %g\n", mat[1][0], mat[1][1], mat[1][2]);
+ // printf("%g %g %g\n", mat[2][0], mat[2][1], mat[2][2]);
+ // printf("%g x^3 + %g x^2 + %g x + %g = 0\n", c[3], c[2], c[1], c[0]);
double imag[3];
FindCubicRoots(c, v, imag);
@@ -522,7 +522,7 @@ void FindCubicRoots(const double coef[4], double real[3], double imag[3])
double d = coef[0];
if(!a || !d){
- // Msg::Error("Degenerate cubic: use a second degree solver!");
+ // Msg::Error("Degenerate cubic: use a second degree solver!");
return;
}
@@ -896,7 +896,8 @@ void signedDistancesPointsTriangle(std::vector<double> &distances,
closePt = p3;
d = sign * std::min(fabs(d), p.distance(p3));
}
- //d = sign * std::min(fabs(d), std::min(std::min(p.distance(p1), p.distance(p2)),p.distance(p3)));
+ //d = sign * std::min(fabs(d), std::min(std::min(p.distance(p1),
+ // p.distance(p2)),p.distance(p3)));
distances[i] = d;
closePts[i] = closePt;
}
@@ -952,41 +953,49 @@ void changeReferential(const int direction,const SPoint3 &p,const SPoint3 &close
const SPoint3 &p1, const SPoint3 &p2, double* xp, double* yp,
double* otherp, double* x, double* y, double* other)
{
- if (direction == 1){
- const SPoint3 &d1=SPoint3(1.0,0.0,0.0);
- const SPoint3 &d=SPoint3(p2.x()-p1.x(),p2.y()-p1.y(),p2.z()-p1.z());
- double norm=sqrt( d.x()*d.x()+d.y()*d.y()+d.z()*d.z() );
- const SPoint3 &dn=SPoint3(d.x()/norm,d.y()/norm,d.z()/norm);
- const SPoint3 &d3=SPoint3(d1.y()*dn.z()-d1.z()*dn.y(),d1.z()*dn.x()-d1.x()*dn.z(),d1.x()*dn.y()-d1.y()*dn.x());
- norm=sqrt( d3.x()*d3.x()+d3.y()*d3.y()+d3.z()*d3.z() );
- const SPoint3 &d3n=SPoint3(d3.x()/norm,d3.y()/norm,d3.z()/norm);
- const SPoint3 &d2=SPoint3(d3n.y()*d1.z()-d3n.z()*d1.y(),d3n.z()*d1.x()-d3n.x()*d1.z(),d3n.x()*d1.y()-d3n.y()*d1.x());
- norm=sqrt( d2.x()*d2.x()+d2.y()*d2.y()+d2.z()*d2.z() );
- const SPoint3 &d2n=SPoint3(d2.x()/norm,d2.y()/norm,d2.z()/norm);
- *xp=p.x()*d1.x()+p.y()*d1.y()+p.z()*d1.z();
- *yp=p.x()*d3n.x()+p.y()*d3n.y()+p.z()*d3n.z();
- *otherp=p.x()*d2n.x()+p.y()*d2n.y()+p.z()*d2n.z();
- *x=closePt.x()*d1.x()+closePt.y()*d1.y()+closePt.z()*d1.z();
- *y=closePt.x()*d3n.x()+closePt.y()*d3n.y()+closePt.z()*d3n.z();
- *other=closePt.x()*d2n.x()+closePt.y()*d2n.y()+closePt.z()*d2n.z();
+ if(direction == 1){
+ const SPoint3 &d1 = SPoint3(1.0, 0.0, 0.0);
+ const SPoint3 &d = SPoint3(p2.x() - p1.x(), p2.y() - p1.y(), p2.z() - p1.z());
+ double norm = sqrt(d.x() * d.x() + d.y() * d.y() + d.z() * d.z());
+ const SPoint3 &dn = SPoint3(d.x() / norm, d.y() / norm, d.z() / norm);
+ const SPoint3 &d3 = SPoint3(d1.y() * dn.z() - d1.z() * dn.y(),
+ d1.z() * dn.x() - d1.x() * dn.z(),
+ d1.x() * dn.y() - d1.y() * dn.x());
+ norm = sqrt(d3.x() * d3.x() + d3.y() * d3.y() + d3.z() * d3.z());
+ const SPoint3 &d3n = SPoint3(d3.x() / norm, d3.y() / norm, d3.z() / norm);
+ const SPoint3 &d2 = SPoint3(d3n.y() * d1.z() - d3n.z() * d1.y(),
+ d3n.z() * d1.x() - d3n.x() * d1.z(),
+ d3n.x() * d1.y() - d3n.y() * d1.x());
+ norm = sqrt(d2.x() * d2.x() + d2.y() * d2.y() + d2.z() * d2.z());
+ const SPoint3 &d2n = SPoint3(d2.x() / norm, d2.y() / norm, d2.z() / norm);
+ *xp = p.x() * d1.x() + p.y() * d1.y() + p.z() * d1.z();
+ *yp = p.x() * d3n.x() + p.y() * d3n.y() + p.z() * d3n.z();
+ *otherp = p.x() * d2n.x() + p.y() * d2n.y() + p.z() * d2n.z();
+ *x = closePt.x() * d1.x() + closePt.y() * d1.y() + closePt.z() * d1.z();
+ *y = closePt.x() * d3n.x() + closePt.y() * d3n.y() + closePt.z() * d3n.z();
+ *other = closePt.x() * d2n.x() + closePt.y() * d2n.y() + closePt.z() * d2n.z();
}
else{
- const SPoint3 &d2=SPoint3(0.0,1.0,0.0);
- const SPoint3 &d=SPoint3(p2.x()-p1.x(),p2.y()-p1.y(),p2.z()-p1.z());
- double norm=sqrt( d.x()*d.x()+d.y()*d.y()+d.z()*d.z() );
- const SPoint3 &dn=SPoint3(d.x()/norm,d.y()/norm,d.z()/norm);
- const SPoint3 &d3=SPoint3(dn.y()*d2.z()-dn.z()*d2.y(),dn.z()*d2.x()-dn.x()*d2.z(),dn.x()*d2.y()-dn.y()*d2.x());
- norm=sqrt( d3.x()*d3.x()+d3.y()*d3.y()+d3.z()*d3.z() );
- const SPoint3 &d3n=SPoint3(d3.x()/norm,d3.y()/norm,d3.z()/norm);
- const SPoint3 &d1=SPoint3(d2.y()*d3n.z()-d2.z()*d3n.y(),d2.z()*d3n.x()-d2.x()*d3n.z(),d2.x()*d3n.y()-d2.y()*d3n.x());
- norm=sqrt( d1.x()*d1.x()+d1.y()*d1.y()+d1.z()*d1.z() );
- const SPoint3 &d1n=SPoint3(d1.x()/norm,d1.y()/norm,d1.z()/norm);
- *xp=p.x()*d2.x()+p.y()*d2.y()+p.z()*d2.z();
- *yp=p.x()*d3n.x()+p.y()*d3n.y()+p.z()*d3n.z();
- *otherp=p.x()*d1n.x()+p.y()*d1n.y()+p.z()*d1n.z();
- *x=closePt.x()*d2.x()+closePt.y()*d2.y()+closePt.z()*d2.z();
- *y=closePt.x()*d3n.x()+closePt.y()*d3n.y()+closePt.z()*d3n.z();
- *other=closePt.x()*d1n.x()+closePt.y()*d1n.y()+closePt.z()*d1n.z();
+ const SPoint3 &d2 = SPoint3(0.0, 1.0, 0.0);
+ const SPoint3 &d = SPoint3(p2.x() - p1.x(), p2.y() - p1.y(), p2.z() - p1.z());
+ double norm = sqrt(d.x() * d.x() + d.y() * d.y() + d.z() * d.z());
+ const SPoint3 &dn = SPoint3(d.x() / norm, d.y() / norm, d.z() / norm);
+ const SPoint3 &d3 = SPoint3(dn.y() * d2.z() - dn.z() * d2.y(),
+ dn.z() * d2.x() - dn.x() * d2.z(),
+ dn.x() * d2.y() - dn.y() * d2.x());
+ norm = sqrt(d3.x() * d3.x() + d3.y() * d3.y() + d3.z() * d3.z());
+ const SPoint3 &d3n = SPoint3(d3.x() / norm, d3.y() / norm, d3.z() / norm);
+ const SPoint3 &d1 = SPoint3(d2.y() * d3n.z() - d2.z() * d3n.y(),
+ d2.z() * d3n.x() - d2.x() * d3n.z(),
+ d2.x() * d3n.y() - d2.y() * d3n.x());
+ norm = sqrt(d1.x() * d1.x() + d1.y() * d1.y() + d1.z() * d1.z());
+ const SPoint3 &d1n = SPoint3(d1.x() / norm, d1.y() / norm, d1.z() / norm);
+ *xp = p.x() * d2.x() + p.y() * d2.y() + p.z() * d2.z();
+ *yp = p.x() * d3n.x() + p.y() * d3n.y() + p.z() * d3n.z();
+ *otherp = p.x() * d1n.x() + p.y() * d1n.y() + p.z() * d1n.z();
+ *x = closePt.x() * d2.x() + closePt.y() * d2.y() + closePt.z() * d2.z();
+ *y = closePt.x() * d3n.x() + closePt.y() * d3n.y() + closePt.z() * d3n.z();
+ *other = closePt.x() * d1n.x() + closePt.y() * d1n.y() + closePt.z() * d1n.z();
}
}
@@ -996,100 +1005,114 @@ int computeDistanceRatio(const double &y, const double &yp, const double &x,
{
double b;
double a;
- if (y==yp){
- b=-y;
- a=0.0;
- }else{
- if (x==xp){
- b=-x;
- a=0.0;
- }else{
- b=(xp*y-x*yp)/(yp-y);
- if (yp==0.0){
+ if (y == yp){
+ b = -y;
+ a = 0.0;
+ }
+ else{
+ if (x == xp){
+ b = -x;
+ a = 0.0;
+ }
+ else{
+ b = (xp * y - x * yp) / (yp - y);
+ if (yp == 0.0){
a=-(b+x)/y;
- }else{
- a=-(b+xp)/yp;
+ }
+ else{
+ a = -(b + xp) / yp;
}
}
}
double ae;
double be;
double ce;
- double da=r1*r1;
- double db=r2*r2;
- if (y==yp){
- ae=1.0/da;
- be=-(2*x)/da;
- ce=(x*x/da)-1.0;
- }else{
- if (x==xp){
- ae = 1.0/db;
- be = -(2.0*y)/db;
- ce = (y*y/db)-1.0;
- }else{
- if (fabs(a)<0.00001){
- ae = 1.0/db;
- be = -(2.0*y)/db;
- ce = (y*y/db)-1.0;
- }else{
- double a2=a*a;
- ae=(1.0/da)+(1.0/(db*a2));
- be=(2.0*y)/(db*a)+(2.0*b)/(a2*db)-((2.0*x)/da);
- ce=(x*x)/da+(b*b)/(db*a2)+(2.0*b*y)/(a*db)+(y*y/db)-1.0;
+ double da = r1 * r1;
+ double db = r2 * r2;
+ if (y == yp){
+ ae = 1.0 / da;
+ be = -(2 * x) / da;
+ ce = (x * x / da) - 1.0;
+ }
+ else{
+ if (x == xp){
+ ae = 1.0 / db;
+ be = -(2.0 * y) / db;
+ ce = (y * y / db) - 1.0;
+ }
+ else{
+ if (fabs(a) < 0.00001){
+ ae = 1.0 / db;
+ be = -(2.0 * y) / db;
+ ce = (y * y / db) - 1.0;
+ }
+ else{
+ double a2 = a * a;
+ ae = (1.0 / da) + (1.0 / (db * a2));
+ be = (2.0 * y)/(db * a) + (2.0 * b) / (a2 * db) - ((2.0 * x) / da);
+ ce = (x * x) / da + (b * b) / (db * a2) +
+ (2.0 * b * y) / (a * db) + (y * y / db) - 1.0;
}
}
}
- double rho=be*be-4*ae*ce;
- double x1,x2,y1,y2,propdist;
- if (rho<0) {
+ double rho = be * be - 4 * ae * ce;
+ double x1, x2, y1, y2, propdist;
+ if (rho < 0) {
return 1;
- }else{
- x1=-(be+sqrt(rho))/(2.0*ae);
- x2=(-be+sqrt(rho))/(2.0*ae);
- if (y==yp){
- y1=-b;
- y2=-b;
- }else{
- if (x==xp){
- y1=x1;
- y2=x2;
- x1=-b;
- x2=-b;
- }else{
- if (fabs(a)<0.00001){
- y1=x1;
- y2=x2;
- x1=-b;
- x2=-b;
- }else{
- y1=-(b+x1)/a;
- y2=-(b+x2)/a;
+ }
+ else{
+ x1 = -(be + sqrt(rho)) / (2.0 * ae);
+ x2 = (-be + sqrt(rho)) / (2.0 * ae);
+ if (y == yp){
+ y1 = -b;
+ y2 = -b;
+ }
+ else{
+ if (x == xp){
+ y1 = x1;
+ y2 = x2;
+ x1 = -b;
+ x2 = -b;
+ }
+ else{
+ if (fabs(a) < 0.00001){
+ y1 = x1;
+ y2 = x2;
+ x1 = -b;
+ x2 = -b;
+ }
+ else{
+ y1 = -(b + x1) / a;
+ y2 = -(b + x2) / a;
}
}
}
- if (x1==x2){
- propdist=(y1-y)/(yp-y);
- if(propdist<0.0){
- propdist=(y2-y)/(yp-y);
+ if (x1 == x2){
+ propdist = (y1 - y) / (yp - y);
+ if(propdist < 0.0){
+ propdist = (y2 - y) / (yp - y);
}
- }else{
- if (xp!=x){
- propdist=(x1-x)/(xp-x);
- if (propdist<0.0){
- propdist=(x2-x)/(xp-x);
+ }
+ else{
+ if (xp != x){
+ propdist = (x1 - x) / (xp - x);
+ if (propdist < 0.0){
+ propdist = (x2 - x) / (xp - x);
}
- }else{
- if (yp!=y){
- propdist=(y1-y)/(yp-y);
- if(propdist<0.0){
- propdist=(y2-y)/(yp-y);
+ }
+ else{
+ if (yp != y){
+ propdist = (y1 - y) / (yp - y);
+ if(propdist < 0.0){
+ propdist = (y2 - y) / (yp - y);
}
- }else{
- propdist=0.01;
+ }
+ else{
+ propdist = 0.01;
}
}
}
- *distance=propdist;
+ *distance = propdist;
return 0;
}
}
@@ -1123,76 +1146,82 @@ void signedDistancesPointsEllipseLine(std::vector<double>&distances,
double xp,yp,x,y,otherp,other,propdist;
if (p1.x()==p2.x()){
direction=1;
- if (fabs(closePt.x()-0.0)<0.00000001) isInYarn[i]=1;
- if (fabs(closePt.x()-2.2)<0.00000001) isInYarn[i]=4;
- if (fabs(closePt.x()-4.4)<0.00000001) isInYarn[i]=2;
- if (fabs(closePt.x()-6.6)<0.00000001) isInYarn[i]=5;
- if (fabs(closePt.x()-8.8)<0.00000001) isInYarn[i]=3;
- if (fabs(closePt.x()-11.0)<0.00000001) isInYarn[i]=1;
- }else{
- if (p1.y()==p2.y()){
- direction=2;
- if (fabs(closePt.y()-0.0)<0.00000001) isInYarn[i]=6;
- if (fabs(closePt.y()-2.2)<0.00000001) isInYarn[i]=7;
- if (fabs(closePt.y()-4.4)<0.00000001) isInYarn[i]=8;
- if (fabs(closePt.y()-6.6)<0.00000001) isInYarn[i]=9;
- if (fabs(closePt.y()-8.8)<0.00000001) isInYarn[i]=10;
- if (fabs(closePt.y()-11.0)<0.00000001) isInYarn[i]=6;
- }else{
- printf("WTF %lf %lf\n",closePt.x(),closePt.y());
+ if (fabs(closePt.x() - 0.0) < 0.00000001) isInYarn[i] = 1;
+ if (fabs(closePt.x() - 2.2) < 0.00000001) isInYarn[i] = 4;
+ if (fabs(closePt.x() - 4.4) < 0.00000001) isInYarn[i] = 2;
+ if (fabs(closePt.x() - 6.6) < 0.00000001) isInYarn[i] = 5;
+ if (fabs(closePt.x() - 8.8) < 0.00000001) isInYarn[i] = 3;
+ if (fabs(closePt.x() - 11.0) < 0.00000001) isInYarn[i] = 1;
+ }
+ else{
+ if (p1.y() == p2.y()){
+ direction = 2;
+ if (fabs(closePt.y() - 0.0) < 0.00000001) isInYarn[i] = 6;
+ if (fabs(closePt.y() - 2.2) < 0.00000001) isInYarn[i] = 7;
+ if (fabs(closePt.y() - 4.4) < 0.00000001) isInYarn[i] = 8;
+ if (fabs(closePt.y() - 6.6) < 0.00000001) isInYarn[i] = 9;
+ if (fabs(closePt.y() - 8.8) < 0.00000001) isInYarn[i] = 10;
+ if (fabs(closePt.y() - 11.0) < 0.00000001) isInYarn[i] = 6;
+ }
+ else{
+ printf("WTF %lf %lf\n", closePt.x(), closePt.y());
}
}
- changeReferential(direction,p,closePt,p1,p2,&xp,&yp,&otherp,&x,&y,&other);
+ changeReferential(direction, p, closePt, p1, p2, &xp, &yp,
+ &otherp, &x, &y, &other);
int result;
- if (fabs(other-otherp)>0.01){
- result=1;
- }else{
- result=computeDistanceRatio(y,yp,x,xp,&propdist,1.1,0.0875);
+ if (fabs(other-otherp) > 0.01){
+ result = 1;
}
- if (result==1){
- distancesE[i] = 1.e10;
- isInYarn[i]=0;
- }else{
- if (propdist<1.0){
- isInYarn[i]=0;
- distancesE[i]=(1.0/propdist)-1.0;
- }else{
- distancesE[i]=(1.0-(1.0/propdist))/3.0;
+ else{
+ result = computeDistanceRatio(y, yp, x, xp, &propdist, 1.1, 0.0875);
+ }
+ if (result == 1){
+ distancesE[i] = 1.e10;
+ isInYarn[i] = 0;
+ }
+ else{
+ if (propdist < 1.0){
+ isInYarn[i] = 0;
+ distancesE[i] = (1.0 / propdist) - 1.0;
+ }
+ else{
+ distancesE[i] = (1.0 - (1.0 / propdist)) / 3.0;
}
}
- }else{
- isInYarn[i]=0;
- distancesE[i]=1000000.0;
+ }
+ else{
+ isInYarn[i] = 0;
+ distancesE[i] = 1000000.0;
}
}
}
-int intersection_segments (SPoint3 &p1, SPoint3 &p2,
- SPoint3 &q1, SPoint3 &q2,
- double x[2])
+int intersection_segments(SPoint3 &p1, SPoint3 &p2,
+ SPoint3 &q1, SPoint3 &q2,
+ double x[2])
{
- double xp_max = std::max(p1.x(),p2.x());
- double yp_max = std::max(p1.y(),p2.y());
- double xq_max = std::max(q1.x(),q2.x());
- double yq_max = std::max(q1.y(),q2.y());
-
- double xp_min = std::min(p1.x(),p2.x());
- double yp_min = std::min(p1.y(),p2.y());
- double xq_min = std::min(q1.x(),q2.x());
- double yq_min = std::min(q1.y(),q2.y());
+ double xp_max = std::max(p1.x(), p2.x());
+ double yp_max = std::max(p1.y(), p2.y());
+ double xq_max = std::max(q1.x(), q2.x());
+ double yq_max = std::max(q1.y(), q2.y());
+
+ double xp_min = std::min(p1.x(), p2.x());
+ double yp_min = std::min(p1.y(), p2.y());
+ double xq_min = std::min(q1.x(), q2.x());
+ double yq_min = std::min(q1.y(), q2.y());
if (yq_min > yp_max || xq_min > xp_max ||
yq_max < yp_min || xq_max < xp_min){
return 0;
}
else{
double A[2][2];
- A[0][0] = p2.x()-p1.x();
- A[0][1] = q1.x()-q2.x();
- A[1][0] = p2.y()-p1.y();
- A[1][1] = q1.y()-q2.y();
- double b[2] = {q1.x()-p1.x(),q1.y()-p1.y()};
- sys2x2(A,b,x);
-
+ A[0][0] = p2.x() - p1.x();
+ A[0][1] = q1.x() - q2.x();
+ A[1][0] = p2.y() - p1.y();
+ A[1][1] = q1.y() - q2.y();
+ double b[2] = {q1.x() - p1.x(), q1.y() - p1.y()};
+ sys2x2(A, b, x);
return (x[0] >= 0.0 && x[0] <= 1. &&
x[1] >= 0.0 && x[1] <= 1.);
}
diff --git a/Plugin/DiscretizationError.cpp b/Plugin/DiscretizationError.cpp
index a4cf71828160b66422c15fce2b6bc1c10d66138f..ca5f416d4f6437efe80cee76ed9a3149eee8bf02 100644
--- a/Plugin/DiscretizationError.cpp
+++ b/Plugin/DiscretizationError.cpp
@@ -26,7 +26,10 @@ extern "C"
std::string GMSH_DiscretizationErrorPlugin::getHelp() const
{
- return "Plugin(DiscretizationError) computes the error between the mesh and the geometry. It does so by supersampling the elements and computing the distance between the supersampled points dans their projection on the geometry.";
+ return "Plugin(DiscretizationError) computes the error between the mesh "
+ "and the geometry. It does so by supersampling the elements and computing "
+ "the distance between the supersampled points dans their projection on "
+ "the geometry.";
}
int GMSH_DiscretizationErrorPlugin::getNbOptions() const
@@ -46,7 +49,8 @@ PView *GMSH_DiscretizationErrorPlugin::execute(PView *v)
double paramQuandt = 1.0 / (nEdgeNodes -1) - 10*tol;
double paramQuandtQuad = 2.0 / (nEdgeNodes - 1) - 10*tol;
int i, j, k, counter;
- double startEstimate[2] = {0.5, 0.5}; // used as a start estimate for u,v when performing an orthogonal projection
+ // used as a start estimate for u,v when performing an orthogonal projection
+ double startEstimate[2] = {0.5, 0.5};
double dx,dy,dz,distance;
std::vector< std::pair<SPoint3,double> > quadDist(nEdgeNodes*nEdgeNodes);
@@ -55,8 +59,8 @@ PView *GMSH_DiscretizationErrorPlugin::execute(PView *v)
PView *v2 = new PView();
PViewDataList *data2 = getDataList(v2);
-
- for (GModel::fiter itFace = GModel::current()->firstFace(); itFace != GModel::current()->lastFace(); ++itFace) {
+ for (GModel::fiter itFace = GModel::current()->firstFace();
+ itFace != GModel::current()->lastFace(); ++itFace) {
// sample quadrangles
/* 13 14 15 16
@@ -64,10 +68,12 @@ PView *GMSH_DiscretizationErrorPlugin::execute(PView *v)
* 5 6 7 8
* 1 2 3 4
*/
- for (std::vector<MQuadrangle*>::iterator itQuad = (*itFace)->quadrangles.begin(); itQuad != (*itFace)->quadrangles.end(); ++itQuad) {
+ for (std::vector<MQuadrangle*>::iterator itQuad = (*itFace)->quadrangles.begin();
+ itQuad != (*itFace)->quadrangles.end(); ++itQuad) {
for (j = 0; j < nEdgeNodes; j++) { // u
for (i = 0; i < nEdgeNodes; i++) { // v
- (*itQuad)->pnt(-1+5*tol+paramQuandtQuad*i, -1+5*tol+paramQuandtQuad*j, 0.0, quadDist[j*(nEdgeNodes) + i].first);
+ (*itQuad)->pnt(-1+5*tol+paramQuandtQuad*i, -1+5*tol+paramQuandtQuad*j,
+ 0.0, quadDist[j*(nEdgeNodes) + i].first);
SPoint3 *point = &quadDist[j*(nEdgeNodes) + i].first;
GPoint closest = (*itFace)->closestPoint(*point,startEstimate);
dx = closest.x() - point->x();
@@ -76,17 +82,18 @@ PView *GMSH_DiscretizationErrorPlugin::execute(PView *v)
quadDist[j*(nEdgeNodes) + i].second = sqrt(dx*dx + dy*dy + dz*dz);
}
}
- for (j = 0; j < nEdgeNodes - 1; j++) {
- for (i = 0; i < nEdgeNodes - 1; i++) {
- std::vector<double> *out = data2->incrementList(1, TYPE_QUA);
- int indices[4] = {i+j*nEdgeNodes, i+j*nEdgeNodes+1, i+(j+1)*nEdgeNodes+1, i+(j+1)*nEdgeNodes};
- for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.x());
- for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.y());
- for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.z());
- for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].second);
+ for (j = 0; j < nEdgeNodes - 1; j++) {
+ for (i = 0; i < nEdgeNodes - 1; i++) {
+ std::vector<double> *out = data2->incrementList(1, TYPE_QUA);
+ int indices[4] = {i+j*nEdgeNodes, i+j*nEdgeNodes+1, i+(j+1)*nEdgeNodes+1,
+ i+(j+1)*nEdgeNodes};
+ for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.x());
+ for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.y());
+ for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].first.z());
+ for (k = 0; k < 4; k++) out->push_back(quadDist[indices[k]].second);
+ }
}
}
- }
// sample triangles
/* 10
@@ -94,11 +101,13 @@ PView *GMSH_DiscretizationErrorPlugin::execute(PView *v)
* 3 5 8
* 1 2 4 7
*/
- for (std::vector<MTriangle*>::iterator itTri = (*itFace)->triangles.begin(); itTri != (*itFace)->triangles.end(); ++itTri) {
+ for (std::vector<MTriangle*>::iterator itTri = (*itFace)->triangles.begin();
+ itTri != (*itFace)->triangles.end(); ++itTri) {
counter = 0;
for (i = 0; i < nEdgeNodes; i++) {
for (j = 0; j < (i + 1); j++) {
- (*itTri)->pnt(5*tol + paramQuandt*(i - j), 5*tol + paramQuandt*j, 0.0, triDist[counter].first);
+ (*itTri)->pnt(5*tol + paramQuandt*(i - j), 5*tol + paramQuandt*j, 0.0,
+ triDist[counter].first);
SPoint3 *point = &triDist[counter].first; // Check : the points are good
GPoint closest = (*itFace)->closestPoint(*point,startEstimate);
dx = (closest.x() - point->x());
@@ -108,39 +117,34 @@ PView *GMSH_DiscretizationErrorPlugin::execute(PView *v)
counter++;
}
}
-
- int indices[3];
- for (j = 0; j < nEdgeNodes - 1; j++) { // row in the triangle
- bool odd = false;
- for (i = 0; i < j*2 + 1; i ++) { // small tri in the row
- if (!odd) {
- indices[0] = i / 2 + (j+1)*j/2;
- indices[1] = i / 2 + (j+1)*j/2 + j+1;
- indices[2] = i / 2 + (j+1)*j/2 + j+2;
- } else {
- indices[0] = (i-1)/2 + (j+1)*j/2;
- indices[1] = (i-1)/2 + (j+1)*j/2 + j+2;
- indices[2] = (i-1)/2 + (j+1)*j/2 + 1;
- }
- std::vector<double> *out = data2->incrementList(1, TYPE_TRI);
- for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.x());
- for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.y());
- for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.z());
- for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].second);
- odd = !odd;
+
+ int indices[3];
+ for (j = 0; j < nEdgeNodes - 1; j++) { // row in the triangle
+ bool odd = false;
+ for (i = 0; i < j*2 + 1; i ++) { // small tri in the row
+ if (!odd) {
+ indices[0] = i / 2 + (j+1)*j/2;
+ indices[1] = i / 2 + (j+1)*j/2 + j+1;
+ indices[2] = i / 2 + (j+1)*j/2 + j+2;
+ }
+ else {
+ indices[0] = (i-1)/2 + (j+1)*j/2;
+ indices[1] = (i-1)/2 + (j+1)*j/2 + j+2;
+ indices[2] = (i-1)/2 + (j+1)*j/2 + 1;
+ }
+ std::vector<double> *out = data2->incrementList(1, TYPE_TRI);
+ for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.x());
+ for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.y());
+ for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].first.z());
+ for (k = 0; k < 3; k++) out->push_back(triDist[indices[k]].second);
+ odd = !odd;
+ }
}
-
- }
-
-
}
-
+
//viusalize stuff
-
-
-
}
-
+
data2->setName("Discretization Error");
data2->setFileName("discretization_err.pos");
data2->finalize();
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index bbac922dc33b638cd3883a175cb3c5f3f93d5d51..8e80559364300f7489991d7d47e779d7321ec9a4 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -20,7 +20,6 @@
#include "laplaceTerm.h"
#include "crossConfTerm.h"
-
StringXNumber DistanceOptions_Number[] = {
{GMSH_FULLRC, "PhysPoint", NULL, 0.},
{GMSH_FULLRC, "PhysLine", NULL, 0.},
@@ -35,7 +34,6 @@ StringXString DistanceOptions_String[] = {
{GMSH_FULLRC, "Filename", NULL, "distance.pos"}
};
-
extern "C"
{
GMSH_Plugin *GMSH_RegisterDistancePlugin()
@@ -43,6 +41,7 @@ extern "C"
return new GMSH_DistancePlugin();
}
}
+
GMSH_DistancePlugin::GMSH_DistancePlugin()
{
_fileName = "text";
@@ -56,14 +55,18 @@ std::string GMSH_DistancePlugin::getHelp() const
{
return "Plugin(Distance) computes distances to physical entities in "
"a mesh.\n\n"
-
- "Define the physical entities to which the distance is computed. If Point=0, Line=0, and Surface=0, then the distance is computed to all the boundaries of the mesh (edges in 2D and faces in 3D).\n\n"
-
- "Computation<0. computes the geometrical euclidian distance (warning: different than the geodesic distance), and Computation=a>0.0 solves a PDE on the mesh with the diffusion constant mu = a*bbox, with bbox being the max size of the bounding box of the mesh (see paper Legrand 2006).\n\n"
-
- "Min Scale and max Scale, scale the distance function. If min Scale<0 and max Scale<0, then no scaling is applied to the distance function.\n\n"
-
- "Plugin(Distance) creates a new distance view and also saves the view in the fileName.pos file.";
+ "Define the physical entities to which the distance is computed. "
+ "If Point=0, Line=0, and Surface=0, then the distance is computed "
+ "to all the boundaries of the mesh (edges in 2D and faces in 3D).\n\n"
+ "Computation<0. computes the geometrical euclidian distance "
+ "(warning: different than the geodesic distance), and Computation=a>0.0 "
+ "solves a PDE on the mesh with the diffusion constant mu = a*bbox, with "
+ "bbox being the max size of the bounding box of the mesh (see paper "
+ "Legrand 2006).\n\n"
+ "Min Scale and max Scale, scale the distance function. If min Scale<0 "
+ "and max Scale<0, then no scaling is applied to the distance function.\n\n"
+ "Plugin(Distance) creates a new distance view and also saves the view "
+ "in the fileName.pos file.";
}
int GMSH_DistancePlugin::getNbOptions() const
@@ -89,7 +92,6 @@ StringXString *GMSH_DistancePlugin::getOptionStr(int iopt)
void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
std::map<MVertex*,double > _distance_map )
{
-
_fileName = DistanceOptions_String[0].def;
_minScale = (double) DistanceOptions_Number[4].def;
_maxScale = (double) DistanceOptions_Number[5].def;
@@ -145,7 +147,6 @@ void GMSH_DistancePlugin::printView(std::vector<GEntity*> _entities,
}
fprintf(fName,"};\n");
fclose(fName);
-
}
PView *GMSH_DistancePlugin::execute(PView *v)
@@ -489,8 +490,6 @@ PView *GMSH_DistancePlugin::execute(PView *v)
#endif
}
-
- //-------------------------------------------------
return view;
}
diff --git a/Post/PViewData.h b/Post/PViewData.h
index c8583262a1d38a6913056ac02dcb37eaa9d373c2..84ad66aae32b8631662d691b670ead57f03a8c8c 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -242,7 +242,8 @@ class PViewData {
// I/O routines
virtual bool writeSTL(std::string fileName);
virtual bool writeTXT(std::string fileName);
- virtual bool writePOS(std::string fileName, bool binary=false, bool parsed=true, bool append=false);
+ virtual bool writePOS(std::string fileName, bool binary=false, bool parsed=true,
+ bool append=false);
virtual bool writeMSH(std::string fileName, bool binary=false, bool savemesh=true);
virtual bool writeMED(std::string fileName);
};
diff --git a/doc/VERSIONS.txt b/doc/VERSIONS.txt
index 2d9b3f60d832530c881e37f7b5923674f1e066ee..a056d129b626e69304154e0f3b615386936b2881 100644
--- a/doc/VERSIONS.txt
+++ b/doc/VERSIONS.txt
@@ -1,12 +1,13 @@
-2.6.0 (): new quadrilateral meshing algoroithms (Blossom and DelQuad);
-new tensor field visualization modes (eigenvectors, ellipsoid, etc.);
-added support for interpolation schemes in .msh file; added support
-for MED3 format; rescale viewport around visible entities (shift+1:1
-in GUI); unified post-processing field export; new experimental
-stereo+camera visualization mode; experimental BAMG & MMG3D support
-for anisotropic mesh generation; new OCC cut&merge algorithm imported
-from Salome; new ability to connect extruded meshes to tetrahedral
-grids using pyramids; Abaqus (INP) mesh export; various bug fixes and
+2.6.0 (): new quadrilateral meshing algorithms (Blossom and DelQuad);
+new solver module based on ONELAB project; new tensor field
+visualization modes (eigenvectors, ellipsoid, etc.); added support for
+interpolation schemes in .msh file; added support for MED3 format;
+rescale viewport around visible entities (shift+1:1 in GUI); unified
+post-processing field export; new experimental stereo+camera
+visualization mode; experimental BAMG & MMG3D support for anisotropic
+mesh generation; new OCC cut&merge algorithm imported from Salome; new
+ability to connect extruded meshes to tetrahedral grids using
+pyramids; Abaqus (INP) mesh export; various bug fixes and
improvements.
2.5.0 (Oct 15, 2010): new compound geometrical entities (for remeshing