diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index e6f021a23723b18ff7387baf81d1cd74697ece55..19591a4342706ac1ea08aaf10ab96633389dc4ab 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -447,7 +447,7 @@ void inline sort2(T &a, T &b)
}
template <class T>
-void sort4(T *t[3])
+void sort4(T *t[4])
{
sort2<T*>(t[0], t[1]);
sort2<T*>(t[2], t[3]);
diff --git a/Mesh/DivideAndConquer.cpp b/Mesh/DivideAndConquer.cpp
index f7f52a6dec8c6363e5a86d4535494a02c3649c08..aa42ddeacd6da401e4ef943f6f6f4caea0801a55 100644
--- a/Mesh/DivideAndConquer.cpp
+++ b/Mesh/DivideAndConquer.cpp
@@ -727,14 +727,14 @@ double DocRecord::Lloyd(int type)
PointRecord &pt = points[i];
std::vector<SPoint2> pts;
voronoiCell (i,pts);
- double E, A;
+ double E = 0., A;
if (!points[i].data){
SPoint2 p (pt.where.h,pt.where.v);
if (type == 0)
- centroidOfPolygon (p,pts, cgs(i,0), cgs(i,1),E, A);
+ centroidOfPolygon(p,pts, cgs(i,0), cgs(i,1), E, A);
else
- centroidOfOrientedBox (pts, 0.0, cgs(i,0),cgs(i,1),E, A);
+ centroidOfOrientedBox(pts, 0.0, cgs(i,0), cgs(i,1), E, A);
}
inertia_tot += E;
}
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index 17f951430c1de952f682ced70bd2ad6c4b72724e..3356af0dc8590f16cd7ad24ef0d50d8fd2c54e27 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -35,10 +35,10 @@ void Frame_field::init_region(GRegion* gr){
std::list<GFace*> faces;
STensor3 m(1.0);
- Nearest_point::init_region(gr);
-
- faces = gr->faces();
-
+ Nearest_point::init_region(gr);
+
+ faces = gr->faces();
+
temp.clear();
field.clear();
@@ -66,7 +66,7 @@ void Frame_field::init_region(GRegion* gr){
void Frame_field::init_face(GFace* gf){
// Fills the auxiliary std::map "temp" with a pair <vertex, STensor3>
- // for each vertex of the face gf.
+ // for each vertex of the face gf.
unsigned int i;
int j;
bool ok;
@@ -83,30 +83,30 @@ void Frame_field::init_face(GFace* gf){
for(i=0;i<gf->getNumMeshElements();i++){
element = gf->getMeshElement(i);
-
+
average_x = 0.0;
average_y = 0.0;
-
+
for(j=0;j<element->getNumVertices();j++){
vertex = element->getVertex(j);
reparamMeshVertexOnFace(vertex,gf,point);
average_x = average_x + point.x();
average_y = average_y + point.y();
- }
-
+ }
+
average_x = average_x/element->getNumVertices();
average_y = average_y/element->getNumVertices();
-
+
for(j=0;j<element->getNumVertices();j++){
vertex = element->getVertex(j);
-
+
if(gf->geomType()==GEntity::CompoundSurface){
ok = translate(gf,octree,vertex,SPoint2(average_x,average_y),v1,v2);
}
else{
ok = improved_translate(gf,vertex,v1,v2);
}
-
+
if(ok){
v1.normalize();
v2.normalize();
@@ -169,7 +169,7 @@ bool Frame_field::translate(GFace* gf,MElementOctree* octree,MVertex* vertex,SPo
}
ok = true; //?
-
+
if(ok){
gp1 = gf->point(x1,y1);
gp2 = gf->point(x2,y2);
@@ -190,27 +190,27 @@ bool Frame_field::improved_translate(GFace* gf,MVertex* vertex,SVector3& v1,SVec
SVector3 s1,s2;
SVector3 normal;
SVector3 basis_u,basis_v;
- Pair<SVector3,SVector3> derivatives;
-
+ Pair<SVector3,SVector3> derivatives;
+
reparamMeshVertexOnFace(vertex,gf,point);
x = point.x();
y = point.y();
-
+
angle = backgroundMesh::current()->getAngle(x,y,0.0);
derivatives = gf->firstDer(point);
-
+
s1 = derivatives.first();
s2 = derivatives.second();
normal = crossprod(s1,s2);
-
+
basis_u = s1;
basis_u.normalize();
basis_v = crossprod(normal,basis_u);
basis_v.normalize();
-
+
v1 = basis_u*cos(angle) + basis_v*sin(angle);
v2 = crossprod(v1,normal);
-
+
return 1;
}
@@ -250,16 +250,16 @@ STensor3 Frame_field::combine(double x,double y,double z){
SVector3 vec1,vec2,vec3;
SVector3 final1,final2;
STensor3 m(1.0),m2(1.0);
-
+
m = search(x,y,z);
m2 = m;
- ok = Nearest_point::search(x,y,z,vec);
-
+ ok = Nearest_point::search(x,y,z,vec);
+
if(ok){
vec1 = SVector3(m.get_m11(),m.get_m21(),m.get_m31());
vec2 = SVector3(m.get_m12(),m.get_m22(),m.get_m32());
vec3 = SVector3(m.get_m13(),m.get_m23(),m.get_m33());
-
+
val1 = fabs(dot(vec,vec1));
val2 = fabs(dot(vec,vec2));
val3 = fabs(dot(vec,vec3));
@@ -273,11 +273,11 @@ STensor3 Frame_field::combine(double x,double y,double z){
else{
other = vec3;
}
-
+
final1 = crossprod(vec,other);
final1.normalize();
final2 = crossprod(vec,final1);
-
+
m2.set_m11(vec.x());
m2.set_m21(vec.y());
m2.set_m31(vec.z());
@@ -288,7 +288,7 @@ STensor3 Frame_field::combine(double x,double y,double z){
m2.set_m23(final2.y());
m2.set_m33(final2.z());
}
-
+
return m2;
}
@@ -303,7 +303,7 @@ double Frame_field::get_ratio(GFace* gf,SPoint2 point){
double val;
double uv[2];
double tab[3];
-
+
uv[0] = point.x();
uv[1] = point.y();
buildMetric(gf,uv,tab);
@@ -427,7 +427,7 @@ GRegion* Frame_field::test(){
}
void Frame_field::clear(){
- Nearest_point::clear();
+ Nearest_point::clear();
temp.clear();
field.clear();
#if defined(HAVE_ANN)
@@ -441,7 +441,7 @@ void Frame_field::clear(){
#endif
}
-// BARYCENTRIC
+// BARYCENTRIC
#include "cross3D.h"
@@ -450,7 +450,7 @@ double min(const double a, const double b) { return (b<a)?b:a; }
double squ(const double a) { return a*a; }
int Frame_field::build_vertex_to_vertices(GEntity* gr, int onWhat, bool initialize){
- // build vertex to vertices data
+ // build vertex to vertices data
std::set<MVertex*> vertices;
if(initialize) vertex_to_vertices.clear();
for(unsigned int i=0; i<gr->getNumMeshElements(); i++){
@@ -471,7 +471,7 @@ int Frame_field::build_vertex_to_vertices(GEntity* gr, int onWhat, bool initiali
vertices.insert(pElem->getVertex((j+k) % n));
vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(pVertex,vertices));
}
-
+
}
}
return vertex_to_vertices.size();
@@ -559,14 +559,14 @@ int Frame_field::findAnnIndex(SPoint3 p){
}
void Frame_field::initFace(GFace* gf){
- // align crosses of gf with the normal (average on neighbour elements)
+ // align crosses of gf with the normal (average on neighbour elements)
// at all vertices of the GFace gf
build_vertex_to_elements(gf);
for(std::map<MVertex*, std::set<MElement*> >::const_iterator it
= vertex_to_elements.begin(); it != vertex_to_elements.end(); it++){
std::set<MElement*> elements = it->second;
SVector3 Area = SVector3(0,0,0);
- for(std::set<MElement*>::const_iterator iter=elements.begin();
+ for(std::set<MElement*>::const_iterator iter=elements.begin();
iter != elements.end(); iter++){
MElement *pElem = *iter;
int n = pElem->getNumVertices();
@@ -621,7 +621,7 @@ void Frame_field::initFace(GFace* gf){
MVertex* pVertex = it->first;
std::set<MElement*> elements = it->second;
if(elements.size() != 2){
- std::cout << "The face has an open boundary" << std::endl;
+ std::cout << "The face has an open boundary" << std::endl;
exit(1);
}
MEdge edg1 = (*elements.begin())->getEdge(0);
@@ -631,7 +631,7 @@ void Frame_field::initFace(GFace* gf){
std::map<MVertex*, STensor3>::iterator iter = crossField.find(pVertex);
if(iter == crossField.end()) {
- std::cout << "This should not happen: cross not found 1" << std::endl;
+ std::cout << "This should not happen: cross not found 1" << std::endl;
exit(1);
}
cross3D x(cross3D((*iter).second));
@@ -650,7 +650,7 @@ void Frame_field::initFace(GFace* gf){
// fixme: we have to rebuild the vertex_to_elements list
// to have a working iterator oin the verticies of gf
// The natural iterator gf->getMeshVertex(i) seems not to work
- // when the option PackingOfParallelogram is on.
+ // when the option PackingOfParallelogram is on.
build_vertex_to_elements(gf);
for(std::map<MVertex*, std::set<MElement*> >::const_iterator it
= vertex_to_elements.begin(); it != vertex_to_elements.end(); it++){
@@ -664,7 +664,7 @@ void Frame_field::initFace(GFace* gf){
cross3D y;
if(iter == crossField.end()){
std::cout << "This should not happen: cross not found 2" << std::endl;
- std::cout << pVertex0->x() << "," << pVertex0->y() << "," << pVertex0->z() << std::endl;
+ std::cout << pVertex0->x() << "," << pVertex0->y() << "," << pVertex0->z() << std::endl;
//exit(1);
y = cross3D();
}
@@ -676,7 +676,7 @@ void Frame_field::initFace(GFace* gf){
MVertex* pVertex = listVertices[index]; //nearest vertex on contour
iter = crossField.find(pVertex);
if(iter == crossField.end()){
- std::cout << "This should not happen: cross not found 3" << std::endl;
+ std::cout << "This should not happen: cross not found 3" << std::endl;
exit(1);
}
cross3D x = cross3D((*iter).second); //nearest cross on contour, x.getFrst() is the normal
@@ -686,7 +686,7 @@ void Frame_field::initFace(GFace* gf){
STensor3 m = convert( y.rotate(conj(x.rotationToOnSurf(y)))); //rotation around the normal
SVector3 v2 = y.getFrst();
if(fabs(angle(v1,v2)) > 1e-8){
- std::cout << "This should not happen: rotation affected the normal" << std::endl;
+ std::cout << "This should not happen: rotation affected the normal" << std::endl;
exit(1);
}
crossField[pVertex0] = m;
@@ -696,7 +696,7 @@ void Frame_field::initFace(GFace* gf){
}
void Frame_field::initRegion(GRegion* gr, int n){
- std::list<GFace*> faces = gr->faces();
+ std::list<GFace*> faces = gr->faces();
for(std::list<GFace*>::const_iterator iter=faces.begin(); iter!=faces.end(); iter++){
initFace(*iter);
// smoothing must be done immediately because crosses on the contour vertices
@@ -734,7 +734,7 @@ double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::cons
MVertex* pVertex0 = iter->first;
std::set<MVertex*> list = iter->second;
cross3D x(m0);
- if(debug) std::cout << "# " << pVertex0->getNum()
+ if(debug) std::cout << "# " << pVertex0->getNum()
<< " with " << list.size() << " neighbours" << std::endl;
SVector3 T = SVector3(0.), dT;
@@ -742,7 +742,7 @@ double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::cons
energy = 0;
for (std::set<MVertex*>::const_iterator it=list.begin(); it != list.end(); ++it){
MVertex *pVertex = *it;
- if(pVertex->getNum() == pVertex0->getNum())
+ if(pVertex->getNum() == pVertex0->getNum())
std::cout << "This should not happen!" << std::endl;
std::map<MVertex*, STensor3>::const_iterator itB = crossField.find(pVertex);
if(itB == crossField.end()){ // not found
@@ -807,7 +807,7 @@ double Frame_field::smooth(){
double enew, eold;
double energy = 0;
- for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
+ for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
= vertex_to_vertices.begin(); iter != vertex_to_vertices.end(); ++iter){
//MVertex* pVertex0 = iter->first;
@@ -832,7 +832,7 @@ double Frame_field::smooth(){
return energy;
}
-void Frame_field::save(const std::vector<std::pair<SPoint3, STensor3> > data,
+void Frame_field::save(const std::vector<std::pair<SPoint3, STensor3> > data,
const std::string& filename){
const cross3D origin(SVector3(1,0,0), SVector3(0,1,0));
SPoint3 p1;
@@ -886,7 +886,7 @@ void Frame_field::recur_connect_vert(MVertex *v,
it != v2v.upper_bound(v) ; ++it){
MVertex *nextV = it->second;
- //change orientation
+ //change orientation
///------------------
printf("*********************** Changing orientation \n");
@@ -900,7 +900,7 @@ void Frame_field::recur_connect_vert(MVertex *v,
nextCross.print("nextCross ");
prod.print("product");
fullMatrix<double> mat(3,3); prod.getMat(mat);
-
+
//find biggest dot product
fullVector<int> Id(3);
for (int i = 0; i < 3; i++){
@@ -913,38 +913,38 @@ void Frame_field::recur_connect_vert(MVertex *v,
}
}
}
-
+
if (Id(0) +Id(1)+ Id(2) != 3 || (Id(0) == 1 && Id(1)==1 && Id(2)==1)) {
- std::cout << "This should not happen: sum should be 0+1+2" << std::endl;
+ std::cout << "This should not happen: sum should be 0+1+2" << std::endl;
printf("Id =%d %d %d \n", Id(0), Id(1), Id(2));
- //exit(1);
+ //exit(1);
return;
}
//create new cross
- fullMatrix<double> newmat(3,3);
+ fullMatrix<double> newmat(3,3);
for (int i = 0; i < 3; i++){
for (int j = 0; j < 3; j++){
newmat(i,j) = nextCross(Id(j),i);
}
}
-
- STensor3 newcross;
+
+ STensor3 newcross;
newcross.setMat(newmat);
crossField[iter->first] = newcross;
- newcross.print("newcross");
+ newcross.print("newcross");
//continue recursion
- recur_connect_vert (nextV, newcross, v2v,touched);
+ recur_connect_vert (nextV, newcross, v2v,touched);
}
}
void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){
-
+
printf("continuous cross field \n");
-
+
//start from a vertex of a face
- std::list<GFace*> faces = gr->faces();
+ std::list<GFace*> faces = gr->faces();
bool foundFace = false;
for(std::list<GFace*>::const_iterator iter=faces.begin(); iter!=faces.end(); iter++){
if (*iter == gf){
@@ -953,14 +953,14 @@ void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){
}
}
if (!foundFace){
- std::cout << "This should not happen: face does not belong to region" << std::endl;
+ std::cout << "This should not happen: face does not belong to region" << std::endl;
exit(1);
}
//build connectivity
build_vertex_to_vertices(gr, -1, true);
std::multimap<MVertex*,MVertex*> v2v;
- for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
+ for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
= vertex_to_vertices.begin(); iter != vertex_to_vertices.end(); ++iter){
MVertex *v = iter->first;
std::set<MVertex*> mySet = iter->second;
@@ -975,7 +975,7 @@ void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){
std::map<MVertex*, STensor3>::iterator iter = crossField.find(beginV);
STensor3 bCross = iter->second;
recur_connect_vert (beginV,bCross,v2v,touched);
-
+
}
void Frame_field::saveCrossField(const std::string& filename, double scale, bool full){
@@ -984,7 +984,7 @@ void Frame_field::saveCrossField(const std::string& filename, double scale, bool
double k = scale;
std::ofstream file(filename.c_str());
file << "View \"cross field\" {\n";
- for(std::map<MVertex*, STensor3>::const_iterator it = crossField.begin();
+ for(std::map<MVertex*, STensor3>::const_iterator it = crossField.begin();
it != crossField.end(); it++){
SPoint3 p = it->first->point();
STensor3 m = it->second;
@@ -992,7 +992,7 @@ void Frame_field::saveCrossField(const std::string& filename, double scale, bool
//double val1=1.0, val2=1.0;
p1 = SPoint3(p.x() + k*m.get_m11(),
p.y() + k*m.get_m21(),
- p.z() + k*m.get_m31());
+ p.z() + k*m.get_m31());
print_segment(p,p1,val1,val2,file);
p1 = SPoint3(p.x() - k*m.get_m11(),
p.y() - k*m.get_m21(),
@@ -1025,7 +1025,7 @@ void Frame_field::save_dist(const std::string& filename){
std::ofstream file(filename.c_str());
file << "View \"Distance\" {\n";
- for(std::map<MEdge, double>::iterator it = crossDist.begin();
+ for(std::map<MEdge, double>::iterator it = crossDist.begin();
it != crossDist.end(); it++){
MVertex* pVerta = it->first.getVertex(0);
MVertex* pVertb = it->first.getVertex(1);
@@ -1076,9 +1076,9 @@ void Frame_field::save_energy(GRegion* gr, const std::string& filename){
//std::vector<double> *out = data->incrementList(1, TYPE_TET, NumNodes);
std::vector<double> *out = data->incrementList(3, TYPE_TET, NumNodes);
for(int j = 0; j < pTet->getNumVertices(); j++)
- out->push_back(pTet->getVertex(j)->x());
+ out->push_back(pTet->getVertex(j)->x());
for(int j = 0; j < pTet->getNumVertices(); j++)
- out->push_back(pTet->getVertex(j)->y());
+ out->push_back(pTet->getVertex(j)->y());
for(int j = 0; j < pTet->getNumVertices(); j++)
out->push_back(pTet->getVertex(j)->z());
for(int j = 0; j < pTet->getNumVertices(); j++){
@@ -1090,7 +1090,7 @@ void Frame_field::save_energy(GRegion* gr, const std::string& filename){
double jac[3][3], inv[3][3];
pTet->getJacobian(u, v, w, jac);
inv3x3(jac, inv);
- SVector3 sum(0,0,0);
+ SVector3 sum(0,0,0);
for(int k = 0; k < pTet->getNumEdges(); k++){
int nod1 = pTet->edges_tetra(k,0);
int nod2 = pTet->edges_tetra(k,1);
@@ -1135,7 +1135,7 @@ void Size_field::init_region(GRegion* gr){
faces = gr->faces();
- boundary.clear();
+ boundary.clear();
for(it=faces.begin();it!=faces.end();it++){
gf = *it;
@@ -1176,7 +1176,7 @@ void Size_field::solve(GRegion* gr){
size_t i;
int count;
- int count2;
+ int count2;
double val;
double volume;
std::set<MVertex*> interior;
@@ -1193,14 +1193,14 @@ void Size_field::solve(GRegion* gr){
count++;
}
//printf("\n");
- //printf("number of boundary vertices = %d\n",count);
+ //printf("number of boundary vertices = %d\n",count);
for(i=0;i<gr->tetrahedra.size();i++){
interior.insert(gr->tetrahedra[i]->getVertex(0));
interior.insert(gr->tetrahedra[i]->getVertex(1));
interior.insert(gr->tetrahedra[i]->getVertex(2));
interior.insert(gr->tetrahedra[i]->getVertex(3));
- }
+ }
for(it=interior.begin();it!=interior.end();it++){
it2 = boundary.find(*it);
@@ -1223,7 +1223,7 @@ void Size_field::solve(GRegion* gr){
count2++;
volume = volume + gr->tetrahedra[i]->getVolume();
}
- //printf("number of tetrahedra = %d\n",count2);
+ //printf("number of tetrahedra = %d\n",count2);
//printf("volume = %f\n",volume);
if(assembler.sizeOfR()){
@@ -1236,7 +1236,7 @@ void Size_field::solve(GRegion* gr){
}
delete system;
- #endif
+ #endif
}
double Size_field::search(double x,double y,double z){
@@ -1249,11 +1249,11 @@ double Size_field::search(double x,double y,double z){
std::map<MVertex*,double>::iterator it2;
std::map<MVertex*,double>::iterator it3;
std::map<MVertex*,double>::iterator it4;
-
- val = 1.0;
-
+
+ val = 1.0;
+
element = (MElement*)octree->find(x,y,z,3,true);
-
+
if(element!=NULL){
temp1[0] = x;
temp1[1] = y;
@@ -1265,7 +1265,7 @@ double Size_field::search(double x,double y,double z){
v = temp2[1];
w = temp2[2];
- it1 = boundary.find(element->getVertex(0));
+ it1 = boundary.find(element->getVertex(0));
it2 = boundary.find(element->getVertex(1));
it3 = boundary.find(element->getVertex(2));
it4 = boundary.find(element->getVertex(3));
@@ -1323,16 +1323,16 @@ void Size_field::print_field(GRegion* gr){
//printf("x = %f, y = %f, z = %f, mesh size = %f\n",x,y,z,it->second);
}
- //printf("\n");
+ //printf("\n");
printf("min mesh size = %f\n",min);
- printf("max mesh size = %f\n",max);
+ printf("max mesh size = %f\n",max);
printf("total number of vertices = %zu\n",boundary.size());
- printf("\n");
+ printf("\n");
std::ofstream file("laplace.pos");
- file << "View \"test\" {\n";
+ file << "View \"test\" {\n";
for(i=0;i<gr->tetrahedra.size();i++){
x1 = gr->tetrahedra[i]->getVertex(0)->x();
@@ -1362,7 +1362,7 @@ void Size_field::print_field(GRegion* gr){
h4 = it4->second;
file << "SS ("
- << x1 << ", " << y1 << ", " << z1 << ", "
+ << x1 << ", " << y1 << ", " << z1 << ", "
<< x2 << ", " << y2 << ", " << z2 << ", "
<< x3 << ", " << y3 << ", " << z3 << ", "
<< x4 << ", " << y4 << ", " << z4 << "){"
@@ -1410,45 +1410,45 @@ void Nearest_point::init_region(GRegion* gr){
std::set<MVertex*>::iterator it2;
fullMatrix<double> gauss_points;
fullVector<double> gauss_weights;
-
+
gaussIntegration::getTriangle(8,gauss_points,gauss_weights);
- gauss_num = gauss_points.size1();
-
+ gauss_num = gauss_points.size1();
+
faces = gr->faces();
-
+
field.clear();
vicinity.clear();
temp.clear();
-
+
for(it=faces.begin();it!=faces.end();it++){
gf = *it;
for(i=0;i<gf->getNumMeshElements();i++){
element = gf->getMeshElement(i);
-
+
x1 = element->getVertex(0)->x();
y1 = element->getVertex(0)->y();
z1 = element->getVertex(0)->z();
-
+
x2 = element->getVertex(1)->x();
y2 = element->getVertex(1)->y();
z2 = element->getVertex(1)->z();
-
+
x3 = element->getVertex(2)->x();
y3 = element->getVertex(2)->y();
z3 = element->getVertex(2)->z();
-
+
for(j=0;j<gauss_num;j++){
u = gauss_points(j,0);
v = gauss_points(j,1);
-
+
x = T(u,v,x1,x2,x3);
y = T(u,v,y1,y2,y3);
z = T(u,v,z1,z2,z3);
-
+
field.push_back(SPoint3(x,y,z));
vicinity.push_back(element);
}
-
+
temp.insert(element->getVertex(0));
temp.insert(element->getVertex(1));
temp.insert(element->getVertex(2));
@@ -1462,9 +1462,9 @@ void Nearest_point::init_region(GRegion* gr){
//field.push_back(SPoint3(x,y,z));
//vicinity.push_back(NULL);
}
-
+
duplicate = annAllocPts(field.size(),3);
-
+
for(i=0;i<field.size();i++){
duplicate[i][0] = field[i].x();
duplicate[i][1] = field[i].y();
@@ -1475,29 +1475,30 @@ void Nearest_point::init_region(GRegion* gr){
#endif
}
-bool Nearest_point::search(double x,double y,double z,SVector3& vec){
+bool Nearest_point::search(double x,double y,double z,SVector3& vec)
+{
int index;
- bool val;
- #if defined(HAVE_ANN)
+ bool val = false;
+#if defined(HAVE_ANN)
double e;
double e2;
SPoint3 found;
ANNpoint query;
ANNidxArray indices;
ANNdistArray distances;
-
+
query = annAllocPt(3);
query[0] = x;
query[1] = y;
query[2] = z;
-
+
indices = new ANNidx[1];
distances = new ANNdist[1];
-
+
e = 0.0;
kd_tree->annkSearch(query,1,indices,distances,e);
index = indices[0];
-
+
annDeallocPt(query);
delete[] indices;
delete[] distances;
@@ -1508,9 +1509,9 @@ bool Nearest_point::search(double x,double y,double z,SVector3& vec){
else{
found = field[index];
}
-
- e2 = 0.000001;
-
+
+ e2 = 0.000001;
+
if(fabs(found.x()-x)>e2 || fabs(found.y()-y)>e2 || fabs(found.z()-z)>e2){
vec = SVector3(found.x()-x,found.y()-y,found.z()-z);
vec.normalize();
@@ -1520,8 +1521,8 @@ bool Nearest_point::search(double x,double y,double z,SVector3& vec){
vec = SVector3(1.0,0.0,0.0);
val = 0;
}
- #endif
-
+#endif
+
return val;
}
@@ -1531,24 +1532,24 @@ double Nearest_point::T(double u,double v,double val1,double val2,double val3){
//The following method comes from this page : gamedev.net/topic/552906-closest-point-on-triangle
//It can also be found on this page : geometrictools.com/LibMathematics/Distance/Distance.html
-SPoint3 Nearest_point::closest(MElement* element,SPoint3 point){
+SPoint3 Nearest_point::closest(MElement* element,SPoint3 point){
SVector3 edge0 = SVector3(element->getVertex(1)->x()-element->getVertex(0)->x(),element->getVertex(1)->y()-element->getVertex(0)->y(),
element->getVertex(1)->z()-element->getVertex(0)->z());
SVector3 edge1 = SVector3(element->getVertex(2)->x()-element->getVertex(0)->x(),element->getVertex(2)->y()-element->getVertex(0)->y(),
element->getVertex(2)->z()-element->getVertex(0)->z());
SVector3 v0 = SVector3(element->getVertex(0)->x()-point.x(),element->getVertex(0)->y()-point.y(),
element->getVertex(0)->z()-point.z());
-
+
double a = dot(edge0,edge0);
double b = dot(edge0,edge1);
double c = dot(edge1,edge1);
double d = dot(edge0,v0);
double e = dot(edge1,v0);
-
+
double det = a*c-b*b;
double s = b*e-c*d;
double t = b*d-a*e;
-
+
if(s+t<det){
if(s<0.0){
if(t<0.0){
@@ -1610,23 +1611,23 @@ SPoint3 Nearest_point::closest(MElement* element,SPoint3 point){
t = 1.0-s;
}
}
-
+
return SPoint3(element->getVertex(0)->x()+s*edge0.x()+t*edge1.x(),element->getVertex(0)->y()+s*edge0.y()+t*edge1.y(),
element->getVertex(0)->z()+s*edge0.z()+t*edge1.z());
}
double Nearest_point::clamp(double x,double min,double max){
double val;
-
+
val = x;
-
+
if(val<min){
val = min;
}
else if(val>max){
val = max;
}
-
+
return val;
}
@@ -1639,11 +1640,11 @@ void Nearest_point::print_field(GRegion* gr){
MElement* element;
MVertex* vertex;
SVector3 vec;
-
+
k = 0.05;
std::ofstream file("nearest.pos");
file << "View \"test\" {\n";
-
+
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
for(j=0;j<element->getNumVertices();j++){
@@ -1658,7 +1659,7 @@ void Nearest_point::print_field(GRegion* gr){
}
}
}
-
+
file << "};\n";
}
@@ -1672,9 +1673,9 @@ void Nearest_point::print_segment(SPoint3 p1,SPoint3 p2,std::ofstream& file){
GRegion* Nearest_point::test(){
GRegion* gr;
GModel* model = GModel::current();
-
+
gr = *(model->firstRegion());
-
+
return gr;
}
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 0f41aea38b87db37b4b781236d9d48ccbabd9d54..5a07a44fff99ce0c30671291bd71fc9e0ac38d60 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -67,7 +67,7 @@ inline double computeEdgeLinearLength_new(BDS_Point *p1, BDS_Point *p2, GFace *
double SCALINGU, double SCALINGV)
{
const int nbSb = 10;
- GPoint GP[nbSb-1];
+ GPoint GP[nbSb];
for (int i = 1; i < nbSb; i++){
double xi = (double)i / nbSb;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index b20d9e2797f22852ec282a308f5a169156db7901..b4a929d09091832abdc7b730fcb2fcb247d4e957 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -2459,9 +2459,9 @@ static std::vector<MVertex*> computeBoundingPoints (const std::vector<MElement*>
std::list<MVertex *> oriented;
{
std::list<MEdge>::iterator itsz = border.begin();
- border.erase (itsz);
oriented.push_back(itsz->getVertex(0));
oriented.push_back(itsz->getVertex(1));
+ border.erase (itsz);
}
while (border.size()){
std::list<MVertex*>::iterator itb = oriented.begin();