diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 1f309a4cecb422ee36da0cd12274a16ce8ee2554..b4881526e2fdb10873e1a34c00960eba60ad77ba 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -779,10 +779,10 @@ static void Mesh3D(GModel *m)
}
if(treat_region_ok && (CTX::instance()->mesh.recombine3DAll ||
gr->meshAttributes.recombine3D)){
- if (CTX::instance()->mesh.optimize){
- optimizeMeshGRegionGmsh opt;
- opt(gr);
- }
+ if (CTX::instance()->mesh.optimize){
+ optimizeMeshGRegionGmsh opt;
+ opt(gr);
+ }
double a = Cpu();
if (CTX::instance()->mesh.recombine3DLevel >= 0){
Recombinator rec;
@@ -797,18 +797,18 @@ static void Mesh3D(GModel *m)
CTX::instance()->mesh.recombine3DConformity);
// 0: no pyramid, 1: single-step, 2: two-steps (conforming),
// true: fill non-conformities with trihedra
- // while(LaplaceSmoothing (gr)){
- // }
+ // while(LaplaceSmoothing (gr)){
+ // }
nb_elements_recombination += post.get_nb_elements();
nb_hexa_recombination += post.get_nb_hexahedra();
vol_element_recombination += post.get_vol_elements();
vol_hexa_recombination += post.get_vol_hexahedra();
// partial export
- // stringstream ss;
- // ss << "yamakawa_part_";
- // ss << gr->tag();
- // ss << ".msh";
- // export_gregion_mesh(gr, ss.str().c_str());
+ // stringstream ss;
+ // ss << "yamakawa_part_";
+ // ss << gr->tag();
+ // ss << ".msh";
+ // export_gregion_mesh(gr, ss.str().c_str());
time_recombination += (Cpu() - a);
}
}
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index ad92e9d78623ef3d0eb8423611522423efb196f0..d35020558aa5b4526ed319b9be2a2b5ead584999 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -43,15 +43,15 @@ void Frame_field::init_region(GRegion* gr){
for(it=faces.begin();it!=faces.end();it++){
gf = *it;
- init_face(gf);
+ init_face(gf);
}
ANNpointArray duplicate = annAllocPts(field.size(),3);
for(i=0;i<field.size();i++){
duplicate[i][0] = field[i].first.x();
- duplicate[i][1] = field[i].first.y();
- duplicate[i][2] = field[i].first.z();
+ duplicate[i][1] = field[i].first.y();
+ duplicate[i][2] = field[i].first.z();
}
kd_tree = new ANNkd_tree(duplicate,field.size(),3);
@@ -70,24 +70,24 @@ void Frame_field::init_face(GFace* gf){
for(i=0;i<gf->storage1.size();i++){
point = gf->storage1[i];
- v1 = gf->storage2[i];
- v2 = gf->storage3[i];
-
- v1.normalize();
- v2.normalize();
- v3 = crossprod(v1,v2);
- v3.normalize();
- m.set_m11(v1.x());
- m.set_m21(v1.y());
- m.set_m31(v1.z());
- m.set_m12(v2.x());
- m.set_m22(v2.y());
- m.set_m32(v2.z());
- m.set_m13(v3.x());
- m.set_m23(v3.y());
- m.set_m33(v3.z());
- field.push_back(std::pair<SPoint3,STensor3>(point,m));
- labels.push_back(gf->tag());
+ v1 = gf->storage2[i];
+ v2 = gf->storage3[i];
+
+ v1.normalize();
+ v2.normalize();
+ v3 = crossprod(v1,v2);
+ v3.normalize();
+ m.set_m11(v1.x());
+ m.set_m21(v1.y());
+ m.set_m31(v1.z());
+ m.set_m12(v2.x());
+ m.set_m22(v2.y());
+ m.set_m32(v2.z());
+ m.set_m13(v3.x());
+ m.set_m23(v3.y());
+ m.set_m33(v3.z());
+ field.push_back(std::pair<SPoint3,STensor3>(point,m));
+ labels.push_back(gf->tag());
}
}
@@ -135,11 +135,11 @@ STensor3 Frame_field::search(double x,double y,double z){
if(fabs(sqrt(distance2)-sqrt(distance1))<e2){
if(labels[index2]<labels[index1]){
- return field[index2].second;
- }
- else{
- return field[index1].second;
- }
+ return field[index2].second;
+ }
+ else{
+ return field[index1].second;
+ }
}
else{
return field[index1].second;
@@ -163,37 +163,37 @@ STensor3 Frame_field::combine(double x,double y,double z){
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));
-
- if(val1<=val2 && val1<=val3){
- other = vec1;
- }
- else if(val2<=val1 && val2<=val3){
- other = vec2;
- }
- else{
- other = vec3;
- }
-
- final1 = crossprod(vec,other);
- final1.normalize();
- final2 = crossprod(vec,final1);
- final2.normalize();
+ 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));
+
+ if(val1<=val2 && val1<=val3){
+ other = vec1;
+ }
+ else if(val2<=val1 && val2<=val3){
+ other = vec2;
+ }
+ else{
+ other = vec3;
+ }
+
+ final1 = crossprod(vec,other);
+ final1.normalize();
+ final2 = crossprod(vec,final1);
+ final2.normalize();
m2.set_m11(vec.x());
- m2.set_m21(vec.y());
- m2.set_m31(vec.z());
- m2.set_m12(final1.x());
- m2.set_m22(final1.y());
- m2.set_m32(final1.z());
- m2.set_m13(final2.x());
- m2.set_m23(final2.y());
- m2.set_m33(final2.z());
+ m2.set_m21(vec.y());
+ m2.set_m31(vec.z());
+ m2.set_m12(final1.x());
+ m2.set_m22(final1.y());
+ m2.set_m32(final1.z());
+ m2.set_m13(final2.x());
+ m2.set_m23(final2.y());
+ m2.set_m33(final2.z());
}
return m2;
@@ -225,33 +225,33 @@ void Frame_field::print_field1(){
for(i=0;i<field.size();i++){
point = field[i].first;
- m = field[i].second;
-
- p1 = SPoint3(point.x() + k*m.get_m11(),
- point.y() + k*m.get_m21(),
- point.z() + k*m.get_m31());
- p2 = SPoint3(point.x() - k*m.get_m11(),
- point.y() - k*m.get_m21(),
- point.z() - k*m.get_m31());
+ m = field[i].second;
+
+ p1 = SPoint3(point.x() + k*m.get_m11(),
+ point.y() + k*m.get_m21(),
+ point.z() + k*m.get_m31());
+ p2 = SPoint3(point.x() - k*m.get_m11(),
+ point.y() - k*m.get_m21(),
+ point.z() - k*m.get_m31());
p3 = SPoint3(point.x() + k*m.get_m12(),
- point.y() + k*m.get_m22(),
- point.z() + k*m.get_m32());
- p4 = SPoint3(point.x() - k*m.get_m12(),
- point.y() - k*m.get_m22(),
- point.z() - k*m.get_m32());
- p5 = SPoint3(point.x() + k*m.get_m13(),
- point.y() + k*m.get_m23(),
- point.z() + k*m.get_m33());
- p6 = SPoint3(point.x() - k*m.get_m13(),
- point.y() - k*m.get_m23(),
- point.z() - k*m.get_m33());
-
- print_segment(point,p1,color1,color2,file);
- print_segment(point,p2,color1,color2,file);
- print_segment(point,p3,color1,color2,file);
- print_segment(point,p4,color1,color2,file);
- print_segment(point,p5,color1,color2,file);
- print_segment(point,p6,color1,color2,file);
+ point.y() + k*m.get_m22(),
+ point.z() + k*m.get_m32());
+ p4 = SPoint3(point.x() - k*m.get_m12(),
+ point.y() - k*m.get_m22(),
+ point.z() - k*m.get_m32());
+ p5 = SPoint3(point.x() + k*m.get_m13(),
+ point.y() + k*m.get_m23(),
+ point.z() + k*m.get_m33());
+ p6 = SPoint3(point.x() - k*m.get_m13(),
+ point.y() - k*m.get_m23(),
+ point.z() - k*m.get_m33());
+
+ print_segment(point,p1,color1,color2,file);
+ print_segment(point,p2,color1,color2,file);
+ print_segment(point,p3,color1,color2,file);
+ print_segment(point,p4,color1,color2,file);
+ print_segment(point,p5,color1,color2,file);
+ print_segment(point,p6,color1,color2,file);
}
file << "};\n";
@@ -278,38 +278,38 @@ void Frame_field::print_field2(GRegion* gr){
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
- for(j=0;j<element->getNumVertices();j++){
- vertex = element->getVertex(j);
- if(vertex->onWhat()->dim()>2){
- point = SPoint3(vertex->x(),vertex->y(),vertex->z());
- m = search(vertex->x(),vertex->y(),vertex->z());
-
- p1 = SPoint3(point.x() + k*m.get_m11(),
- point.y() + k*m.get_m21(),
- point.z() + k*m.get_m31());
- p2 = SPoint3(point.x() - k*m.get_m11(),
- point.y() - k*m.get_m21(),
- point.z() - k*m.get_m31());
- p3 = SPoint3(point.x() + k*m.get_m12(),
- point.y() + k*m.get_m22(),
- point.z() + k*m.get_m32());
- p4 = SPoint3(point.x() - k*m.get_m12(),
- point.y() - k*m.get_m22(),
- point.z() - k*m.get_m32());
- p5 = SPoint3(point.x() + k*m.get_m13(),
- point.y() + k*m.get_m23(),
- point.z() + k*m.get_m33());
- p6 = SPoint3(point.x() - k*m.get_m13(),
- point.y() - k*m.get_m23(),
- point.z() - k*m.get_m33());
-
- print_segment(point,p1,color1,color2,file);
- print_segment(point,p2,color1,color2,file);
- print_segment(point,p3,color1,color2,file);
- print_segment(point,p4,color1,color2,file);
- print_segment(point,p5,color1,color2,file);
- print_segment(point,p6,color1,color2,file);
- }
+ for(j=0;j<element->getNumVertices();j++){
+ vertex = element->getVertex(j);
+ if(vertex->onWhat()->dim()>2){
+ point = SPoint3(vertex->x(),vertex->y(),vertex->z());
+ m = search(vertex->x(),vertex->y(),vertex->z());
+
+ p1 = SPoint3(point.x() + k*m.get_m11(),
+ point.y() + k*m.get_m21(),
+ point.z() + k*m.get_m31());
+ p2 = SPoint3(point.x() - k*m.get_m11(),
+ point.y() - k*m.get_m21(),
+ point.z() - k*m.get_m31());
+ p3 = SPoint3(point.x() + k*m.get_m12(),
+ point.y() + k*m.get_m22(),
+ point.z() + k*m.get_m32());
+ p4 = SPoint3(point.x() - k*m.get_m12(),
+ point.y() - k*m.get_m22(),
+ point.z() - k*m.get_m32());
+ p5 = SPoint3(point.x() + k*m.get_m13(),
+ point.y() + k*m.get_m23(),
+ point.z() + k*m.get_m33());
+ p6 = SPoint3(point.x() - k*m.get_m13(),
+ point.y() - k*m.get_m23(),
+ point.z() - k*m.get_m33());
+
+ print_segment(point,p1,color1,color2,file);
+ print_segment(point,p2,color1,color2,file);
+ print_segment(point,p3,color1,color2,file);
+ print_segment(point,p4,color1,color2,file);
+ print_segment(point,p5,color1,color2,file);
+ print_segment(point,p6,color1,color2,file);
+ }
}
}
@@ -359,14 +359,14 @@ int Frame_field::build_vertex_to_vertices(GEntity* gr, int onWhat, bool initiali
std::map<MVertex*,std::set<MVertex*> >::iterator it = vertex_to_vertices.find(pVertex);
if(it != vertex_to_vertices.end()){
- for(unsigned int k=1; k<n; k++)
- it->second.insert(pElem->getVertex((j+k) % n));
+ for(unsigned int k=1; k<n; k++)
+ it->second.insert(pElem->getVertex((j+k) % n));
}
else{
- vertices.clear();
- for(unsigned int k=1; k<n; k++)
- vertices.insert(pElem->getVertex((j+k) % n));
- vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(pVertex,vertices));
+ vertices.clear();
+ for(unsigned int k=1; k<n; k++)
+ vertices.insert(pElem->getVertex((j+k) % n));
+ vertex_to_vertices.insert(std::pair<MVertex*,std::set<MVertex*> >(pVertex,vertices));
}
}
@@ -384,11 +384,11 @@ int Frame_field::build_vertex_to_elements(GEntity* gr, bool initialize){
MVertex* pVertex = pElem->getVertex(j);
std::map<MVertex*,std::set<MElement*> >::iterator it = vertex_to_elements.find(pVertex);
if(it != vertex_to_elements.end())
- it->second.insert(pElem);
+ it->second.insert(pElem);
else{
- elements.clear();
- elements.insert(pElem);
- vertex_to_elements.insert(std::pair<MVertex*,std::set<MElement*> >(pVertex,elements));
+ elements.clear();
+ elements.insert(pElem);
+ vertex_to_elements.insert(std::pair<MVertex*,std::set<MElement*> >(pVertex,elements));
}
}
}
@@ -402,7 +402,7 @@ void Frame_field::build_listVertices(GEntity* gr, int dim, bool initialize){
for(int j=0; j<pElem->getNumVertices(); j++){
MVertex * pVertex = pElem->getVertex(j);
if(pVertex->onWhat()->dim() == dim)
- list.insert(pVertex);
+ list.insert(pVertex);
}
}
if(initialize) listVertices.clear();
@@ -459,19 +459,19 @@ void Frame_field::initFace(GFace* gf){
// 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++){
+ = 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();
- iter != elements.end(); iter++){
+ iter != elements.end(); iter++){
MElement *pElem = *iter;
int n = pElem->getNumVertices();
int i;
for(i=0; i<n; i++){
- if(pElem->getVertex(i) == it->first) break;
- if(i == n-1) {
- std::cout << "This should not happen" << std:: endl; exit(1);
- }
+ if(pElem->getVertex(i) == it->first) break;
+ if(i == n-1) {
+ std::cout << "This should not happen" << std:: endl; exit(1);
+ }
}
SVector3 V0 = pElem->getVertex(i)->point();
SVector3 V1 = pElem->getVertex((i+n-1)%n)->point(); // previous vertex
@@ -513,7 +513,7 @@ void Frame_field::initFace(GFace* gf){
// align crosses of the contour of "gf" with the tangent to the contour
for(std::map<MVertex*, std::set<MElement*> >::const_iterator it
- = vertex_to_elements.begin(); it != vertex_to_elements.end(); it++){
+ = vertex_to_elements.begin(); it != vertex_to_elements.end(); it++){
MVertex* pVertex = it->first;
std::set<MElement*> elements = it->second;
if(elements.size() != 2){
@@ -549,7 +549,7 @@ void Frame_field::initFace(GFace* gf){
// 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++){
+ = vertex_to_elements.begin(); it != vertex_to_elements.end(); it++){
//for(unsigned int i=0; i<gf->getNumMeshVertices(); i++){
//MVertex* pVertex0 = gf->getMeshVertex(i);
@@ -631,7 +631,7 @@ double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::cons
std::set<MVertex*> list = iter->second;
cross3D x(m0);
if(debug) std::cout << "# " << pVertex0->getNum()
- << " with " << list.size() << " neighbours" << std::endl;
+ << " with " << list.size() << " neighbours" << std::endl;
SVector3 T = SVector3(0.), dT;
double temp = 0.;
@@ -675,87 +675,87 @@ double Frame_field::findBarycenter(std::map<MVertex*, std::set<MVertex*> >::cons
}
void Frame_field::buildSmoothness(){
- GModel *m = GModel::current();
- std::vector<GEntity*> entities;
- m->getEntities(entities);
- //pour commencer on va creer une map de connectique des Mvertex
- std::map<MVertex*, std::vector<MVertex* > > Neighbours;
- for(unsigned int i = 0; i < entities.size(); i++){
- GEntity* eTmp = entities[i];
- for (unsigned int j = 0; j < eTmp->getNumMeshElements();j++){
- MElement* elem = eTmp->getMeshElement(j);
- for (int k = 0;k < elem->getNumVertices();k++){
- for (int l = k;l < elem->getNumVertices();l++){
- if (k != l){
- MVertex* v1 = elem->getVertex(k);
- MVertex* v2 = elem->getVertex(l);
- Neighbours[v1].push_back(v2);
- Neighbours[v2].push_back(v1);
- }
- }
- }
- }
- }
- for(unsigned int i = 0; i < entities.size(); i++){
- for (unsigned int j = 0;j < entities[i]->mesh_vertices.size();j++){
- //on va traiter chaque noeud
- std::set<MVertex*> V1;
- std::set<MVertex*> V2;
- std::set<MVertex*> V3;
- MVertex* v0 = entities[i]->mesh_vertices[j];
- V1.insert(v0);
- std::vector<MVertex* > v0vec = Neighbours[v0];
- for (unsigned int k = 0;k < v0vec.size();k++){
- V1.insert(v0vec[k]);
- }
- for (std::set<MVertex*>::iterator itSet = V1.begin();itSet != V1.end();itSet++){
- MVertex* vTmp = (*itSet);
- V2.insert(vTmp);
- v0vec = Neighbours[vTmp];
- for (unsigned int k = 0;k < v0vec.size();k++){
- V2.insert(v0vec[k]);
- }
- }
- for (std::set<MVertex*>::iterator itSet = V2.begin();itSet != V2.end();itSet++){
- MVertex* vTmp = (*itSet);
- V3.insert(vTmp);
- v0vec = Neighbours[vTmp];
- for (unsigned int k = 0;k < v0vec.size();k++){
- V3.insert(v0vec[k]);
- }
- }
- //we have all three set here, time to compute the smoothnesses for each one
- std::vector<cross3D> C1;
- std::vector<cross3D> C2;
- std::vector<cross3D> C3;
- double S1 = 0.0;
- double S2 = 0.0;
- double S3 = 0.0;
- for (std::set<MVertex*>::iterator itSet = V1.begin();itSet != V1.end();itSet++){
- MVertex* vTmp = (*itSet);
- STensor3 tTmp = crossField[vTmp];
- cross3D cTmp = cross3D(tTmp);
- C1.push_back(cTmp);
- }
- for (std::set<MVertex*>::iterator itSet = V2.begin();itSet != V2.end();itSet++){
- MVertex* vTmp = (*itSet);
- STensor3 tTmp = crossField[vTmp];
- cross3D cTmp = cross3D(tTmp);
- C2.push_back(cTmp);
- }
- for (std::set<MVertex*>::iterator itSet = V3.begin();itSet != V3.end();itSet++){
- MVertex* vTmp = (*itSet);
- STensor3 tTmp = crossField[vTmp];
- cross3D cTmp = cross3D(tTmp);
- C3.push_back(cTmp);
- }
- S1 = computeSetSmoothness(C1);
- S2 = computeSetSmoothness(C2);
- S3 = computeSetSmoothness(C3);
- double finalSmoothness = (9.0 * S1 + 3.0 * S2 + S3) / 13.0 ;
- crossFieldSmoothness[v0] = finalSmoothness;
- }
- }
+ GModel *m = GModel::current();
+ std::vector<GEntity*> entities;
+ m->getEntities(entities);
+ //pour commencer on va creer une map de connectique des Mvertex
+ std::map<MVertex*, std::vector<MVertex* > > Neighbours;
+ for(unsigned int i = 0; i < entities.size(); i++){
+ GEntity* eTmp = entities[i];
+ for (unsigned int j = 0; j < eTmp->getNumMeshElements();j++){
+ MElement* elem = eTmp->getMeshElement(j);
+ for (int k = 0;k < elem->getNumVertices();k++){
+ for (int l = k;l < elem->getNumVertices();l++){
+ if (k != l){
+ MVertex* v1 = elem->getVertex(k);
+ MVertex* v2 = elem->getVertex(l);
+ Neighbours[v1].push_back(v2);
+ Neighbours[v2].push_back(v1);
+ }
+ }
+ }
+ }
+ }
+ for(unsigned int i = 0; i < entities.size(); i++){
+ for (unsigned int j = 0;j < entities[i]->mesh_vertices.size();j++){
+ //on va traiter chaque noeud
+ std::set<MVertex*> V1;
+ std::set<MVertex*> V2;
+ std::set<MVertex*> V3;
+ MVertex* v0 = entities[i]->mesh_vertices[j];
+ V1.insert(v0);
+ std::vector<MVertex* > v0vec = Neighbours[v0];
+ for (unsigned int k = 0;k < v0vec.size();k++){
+ V1.insert(v0vec[k]);
+ }
+ for (std::set<MVertex*>::iterator itSet = V1.begin();itSet != V1.end();itSet++){
+ MVertex* vTmp = (*itSet);
+ V2.insert(vTmp);
+ v0vec = Neighbours[vTmp];
+ for (unsigned int k = 0;k < v0vec.size();k++){
+ V2.insert(v0vec[k]);
+ }
+ }
+ for (std::set<MVertex*>::iterator itSet = V2.begin();itSet != V2.end();itSet++){
+ MVertex* vTmp = (*itSet);
+ V3.insert(vTmp);
+ v0vec = Neighbours[vTmp];
+ for (unsigned int k = 0;k < v0vec.size();k++){
+ V3.insert(v0vec[k]);
+ }
+ }
+ //we have all three set here, time to compute the smoothnesses for each one
+ std::vector<cross3D> C1;
+ std::vector<cross3D> C2;
+ std::vector<cross3D> C3;
+ double S1 = 0.0;
+ double S2 = 0.0;
+ double S3 = 0.0;
+ for (std::set<MVertex*>::iterator itSet = V1.begin();itSet != V1.end();itSet++){
+ MVertex* vTmp = (*itSet);
+ STensor3 tTmp = crossField[vTmp];
+ cross3D cTmp = cross3D(tTmp);
+ C1.push_back(cTmp);
+ }
+ for (std::set<MVertex*>::iterator itSet = V2.begin();itSet != V2.end();itSet++){
+ MVertex* vTmp = (*itSet);
+ STensor3 tTmp = crossField[vTmp];
+ cross3D cTmp = cross3D(tTmp);
+ C2.push_back(cTmp);
+ }
+ for (std::set<MVertex*>::iterator itSet = V3.begin();itSet != V3.end();itSet++){
+ MVertex* vTmp = (*itSet);
+ STensor3 tTmp = crossField[vTmp];
+ cross3D cTmp = cross3D(tTmp);
+ C3.push_back(cTmp);
+ }
+ S1 = computeSetSmoothness(C1);
+ S2 = computeSetSmoothness(C2);
+ S3 = computeSetSmoothness(C3);
+ double finalSmoothness = (9.0 * S1 + 3.0 * S2 + S3) / 13.0 ;
+ crossFieldSmoothness[v0] = finalSmoothness;
+ }
+ }
}
double Frame_field::smoothFace(GFace *gf, int n){
@@ -787,7 +787,7 @@ double Frame_field::smooth(){
double energy = 0;
for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
- = vertex_to_vertices.begin(); iter != vertex_to_vertices.end(); ++iter){
+ = vertex_to_vertices.begin(); iter != vertex_to_vertices.end(); ++iter){
//MVertex* pVertex0 = iter->first;
SVector3 T(0, 0, 0);
@@ -812,7 +812,7 @@ double Frame_field::smooth(){
}
void Frame_field::save(const std::vector<std::pair<SPoint3, STensor3> > data,
- const std::string& filename){
+ const std::string& filename){
const cross3D origin(SVector3(1,0,0), SVector3(0,1,0));
SPoint3 p1;
double k = 0.1;
@@ -823,28 +823,28 @@ void Frame_field::save(const std::vector<std::pair<SPoint3, STensor3> > data,
STensor3 m = data[i].second;
double val1 = eulerAngleFromQtn(cross3D(m).rotationTo(origin)), val2=val1;
p1 = SPoint3(p.x() + k*m.get_m11(),
- p.y() + k*m.get_m21(),
- p.z() + k*m.get_m31());
+ p.y() + k*m.get_m21(),
+ 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(),
- p.z() - k*m.get_m31());
+ p.y() - k*m.get_m21(),
+ p.z() - k*m.get_m31());
print_segment(p,p1,val1,val2,file);
p1 = SPoint3(p.x() + k*m.get_m12(),
- p.y() + k*m.get_m22(),
- p.z() + k*m.get_m32());
+ p.y() + k*m.get_m22(),
+ p.z() + k*m.get_m32());
print_segment(p,p1,val1,val2,file);
p1 = SPoint3(p.x() - k*m.get_m12(),
- p.y() - k*m.get_m22(),
- p.z() - k*m.get_m32());
+ p.y() - k*m.get_m22(),
+ p.z() - k*m.get_m32());
print_segment(p,p1,val1,val2,file);
p1 = SPoint3(p.x() + k*m.get_m13(),
- p.y() + k*m.get_m23(),
- p.z() + k*m.get_m33());
+ p.y() + k*m.get_m23(),
+ p.z() + k*m.get_m33());
print_segment(p,p1,val1,val2,file);
p1 = SPoint3(p.x() - k*m.get_m13(),
- p.y() - k*m.get_m23(),
- p.z() - k*m.get_m33());
+ p.y() - k*m.get_m23(),
+ p.z() - k*m.get_m33());
print_segment(p,p1,val1,val2,file);
}
file << "};\n";
@@ -854,10 +854,10 @@ void Frame_field::save(const std::vector<std::pair<SPoint3, STensor3> > data,
// BARYCENTRIC End
void Frame_field::recur_connect_vert(FILE *fi, int count,
- MVertex *v,
- STensor3 &cross,
- std::multimap<MVertex*,MVertex*> &v2v,
- std::set<MVertex*> &touched)
+ MVertex *v,
+ STensor3 &cross,
+ std::multimap<MVertex*,MVertex*> &v2v,
+ std::set<MVertex*> &touched)
{
if (touched.find(v) != touched.end()) return;
touched.insert(v);
@@ -955,7 +955,7 @@ void Frame_field::continuousCrossField(GRegion *gr, GFace *gf)
build_vertex_to_vertices(gr, -1, true);
std::multimap<MVertex*,MVertex*> v2v;
for(std::map<MVertex*, std::set<MVertex*> >::const_iterator iter
- = vertex_to_vertices.begin(); iter != vertex_to_vertices.end(); ++iter){
+ = vertex_to_vertices.begin(); iter != vertex_to_vertices.end(); ++iter){
MVertex *v = iter->first;
std::set<MVertex*> mySet = iter->second;
for (std::set<MVertex*>::iterator it = mySet.begin(); it!=mySet.end(); ++it){
@@ -996,30 +996,30 @@ void Frame_field::saveCrossField(const std::string& filename, double scale, bool
//double val1 = eulerAngleFromQtn(cross3D(m).rotationTo(origin))*180./M_PI, val2=val1;
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.y() + k*m.get_m21(),
+ 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(),
- p.z() - k*m.get_m31());
+ p.y() - k*m.get_m21(),
+ p.z() - k*m.get_m31());
if(full) print_segment(p,p1,val1,val2,file);
val1=2.0; val2=2.0;
p1 = SPoint3(p.x() + k*m.get_m12(),
- p.y() + k*m.get_m22(),
- p.z() + k*m.get_m32());
+ p.y() + k*m.get_m22(),
+ p.z() + k*m.get_m32());
print_segment(p,p1,val1,val2,file);
p1 = SPoint3(p.x() - k*m.get_m12(),
- p.y() - k*m.get_m22(),
- p.z() - k*m.get_m32());
+ p.y() - k*m.get_m22(),
+ p.z() - k*m.get_m32());
if(full) print_segment(p,p1,val1,val2,file);
val1=3.0; val2=3.0;
p1 = SPoint3(p.x() + k*m.get_m13(),
- p.y() + k*m.get_m23(),
- p.z() + k*m.get_m33());
+ p.y() + k*m.get_m23(),
+ p.z() + k*m.get_m33());
if(full) print_segment(p,p1,val1,val2,file);
p1 = SPoint3(p.x() - k*m.get_m13(),
- p.y() - k*m.get_m23(),
- p.z() - k*m.get_m33());
+ p.y() - k*m.get_m23(),
+ p.z() - k*m.get_m33());
if(full) print_segment(p,p1,val1,val2,file);
}
file << "};\n";
@@ -1038,9 +1038,9 @@ void Frame_field::save_dist(const std::string& filename){
if(it->first.length())
value /= it->first.length();
file << "SL ("
- << pVerta->x() << ", " << pVerta->y() << ", " << pVerta->z() << ", "
- << pVertb->x() << ", " << pVertb->y() << ", " << pVertb->z() << ")"
- << "{" << value << "," << value << "};\n";
+ << pVerta->x() << ", " << pVerta->y() << ", " << pVerta->z() << ", "
+ << pVertb->x() << ", " << pVertb->y() << ", " << pVertb->z() << ")"
+ << "{" << value << "," << value << "};\n";
}
file << "};\n";
file.close();
@@ -1055,9 +1055,9 @@ void Frame_field::checkAnnData(GEntity* ge, const std::string& filename){
MVertex* pVertb = listVertices[findAnnIndex(pVerta->point())];
double value = pVerta->distance(pVertb);
file << "SL ("
- << pVerta->x() << ", " << pVerta->y() << ", " << pVerta->z() << ", "
- << pVertb->x() << ", " << pVertb->y() << ", " << pVertb->z() << ")"
- << "{" << value << "," << value << "};\n";
+ << pVerta->x() << ", " << pVerta->y() << ", " << pVerta->z() << ", "
+ << pVertb->x() << ", " << pVertb->y() << ", " << pVertb->z() << ")"
+ << "{" << value << "," << value << "};\n";
}
file << "};\n";
file.close();
@@ -1076,9 +1076,9 @@ void Frame_field::save_energy(GRegion* gr, const std::string& filename){
for(unsigned int i = 0; i < gr->getNumMeshElements(); i++){
pElem = gr->getMeshElement(i);
MTetrahedron* pTet = new MTetrahedron(pElem->getVertex(0),
- pElem->getVertex(1),
- pElem->getVertex(2),
- pElem->getVertex(3));
+ pElem->getVertex(1),
+ pElem->getVertex(2),
+ pElem->getVertex(3));
//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++)
@@ -1098,15 +1098,15 @@ void Frame_field::save_energy(GRegion* gr, const std::string& filename){
inv3x3(jac, inv);
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);
- double grd1[3], grd2[3];
- matvec(inv, gsf[nod1], grd1);
- matvec(inv, gsf[nod2], grd2);
- SVector3 esf = sf[nod1] * SVector3(grd2) - sf[nod2] * SVector3(grd1);
- std::map<MEdge, double, Less_Edge>::iterator it = crossDist.find(pTet->getEdge(k));
- sum += it->second * esf;
- //sum += (pTet->getVertex(nod2)->z() - pTet->getVertex(nod1)->z()) * esf;
+ int nod1 = pTet->edges_tetra(k,0);
+ int nod2 = pTet->edges_tetra(k,1);
+ double grd1[3], grd2[3];
+ matvec(inv, gsf[nod1], grd1);
+ matvec(inv, gsf[nod2], grd2);
+ SVector3 esf = sf[nod1] * SVector3(grd2) - sf[nod2] * SVector3(grd1);
+ std::map<MEdge, double, Less_Edge>::iterator it = crossDist.find(pTet->getEdge(k));
+ sum += it->second * esf;
+ //sum += (pTet->getVertex(nod2)->z() - pTet->getVertex(nod1)->z()) * esf;
}
out->push_back(sum[0]);
out->push_back(sum[1]);
@@ -1152,20 +1152,20 @@ void Size_field::init_region(GRegion* gr){
for(it=faces.begin();it!=faces.end();it++){
gf = *it;
- for(i=0;i<gf->storage1.size();i++){
- point = gf->storage1[i];
- h = gf->storage4[i];
+ for(i=0;i<gf->storage1.size();i++){
+ point = gf->storage1[i];
+ h = gf->storage4[i];
- field.push_back(std::pair<SPoint3,double>(point,h));
- }
+ field.push_back(std::pair<SPoint3,double>(point,h));
+ }
}
ANNpointArray duplicate = annAllocPts(field.size(),3);
for(i=0;i<field.size();i++){
duplicate[i][0] = field[i].first.x();
- duplicate[i][1] = field[i].first.y();
- duplicate[i][2] = field[i].first.z();
+ duplicate[i][1] = field[i].first.y();
+ duplicate[i][2] = field[i].first.z();
}
kd_tree = new ANNkd_tree(duplicate,field.size(),3);
@@ -1182,22 +1182,22 @@ void Size_field::init_region(GRegion* gr){
for(it=faces.begin();it!=faces.end();it++){
gf = *it;
- for(i=0;i<gf->getNumMeshElements();i++){
+ for(i=0;i<gf->getNumMeshElements();i++){
element = gf->getMeshElement(i);
- for(j=0;j<element->getNumVertices();j++){
- vertex = element->getVertex(j);
+ for(j=0;j<element->getNumVertices();j++){
+ vertex = element->getVertex(j);
- query[0] = vertex->x();
- query[1] = vertex->y();
- query[2] = vertex->z();
+ query[0] = vertex->x();
+ query[1] = vertex->y();
+ query[2] = vertex->z();
- kd_tree->annkSearch(query,1,indices,distances,e);
- index = indices[0];
+ kd_tree->annkSearch(query,1,indices,distances,e);
+ index = indices[0];
- boundary.insert(std::pair<MVertex*,double>(vertex,field[index].second));
- }
- }
+ boundary.insert(std::pair<MVertex*,double>(vertex,field[index].second));
+ }
+ }
}
octree = new MElementOctree(model);
@@ -1456,53 +1456,53 @@ void Nearest_point::init_region(GRegion* gr){
for(it=faces.begin();it!=faces.end();it++){
gf = *it;
- for(i=0;i<gf->getNumMeshElements();i++){
- element = gf->getMeshElement(i);
+ 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();
+ 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();
+ 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();
+ 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);
+ 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);
+ 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);
- }
+ 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));
- }
+ temp.insert(element->getVertex(0));
+ temp.insert(element->getVertex(1));
+ temp.insert(element->getVertex(2));
+ }
}
for(it2=temp.begin();it2!=temp.end();it2++){
- x = (*it2)->x();
- y = (*it2)->y();
- z = (*it2)->z();
+ x = (*it2)->x();
+ y = (*it2)->y();
+ z = (*it2)->z();
//field.push_back(SPoint3(x,y,z));
- //vicinity.push_back(NULL);
+ //vicinity.push_back(NULL);
}
ANNpointArray duplicate = annAllocPts(field.size(),3);
for(i=0;i<field.size();i++){
- duplicate[i][0] = field[i].x();
- duplicate[i][1] = field[i].y();
- duplicate[i][2] = field[i].z();
+ duplicate[i][0] = field[i].x();
+ duplicate[i][1] = field[i].y();
+ duplicate[i][2] = field[i].z();
}
kd_tree = new ANNkd_tree(duplicate,field.size(),3);
@@ -1538,7 +1538,7 @@ bool Nearest_point::search(double x,double y,double z,SVector3& vec)
delete[] distances;
if(vicinity[index]!=NULL){
- found = closest(vicinity[index],SPoint3(x,y,z));
+ found = closest(vicinity[index],SPoint3(x,y,z));
}
else{
found = field[index];
@@ -1549,11 +1549,11 @@ bool Nearest_point::search(double x,double y,double z,SVector3& vec)
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();
- val = 1;
+ val = 1;
}
else{
vec = SVector3(1.0,0.0,0.0);
- val = 0;
+ val = 0;
}
#endif
@@ -1568,11 +1568,11 @@ double Nearest_point::T(double u,double v,double val1,double val2,double val3){
//It can also be found on this page : geometrictools.com/LibMathematics/Distance/Distance.html
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());
+ 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());
+ 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());
+ element->getVertex(0)->z()-point.z());
double a = dot(edge0,edge0);
double b = dot(edge0,edge1);
@@ -1586,68 +1586,68 @@ SPoint3 Nearest_point::closest(MElement* element,SPoint3 point){
if(s+t<det){
if(s<0.0){
- if(t<0.0){
- if(d<0.0){
- s = clamp(-d/a,0.0,1.0);
- t = 0.0;
- }
- else{
- s = 0.0;
- t = clamp(-e/c,0.0,1.0);
- }
- }
- else{
- s = 0.0;
- t = clamp(-e/c,0.0,1.0);
- }
- }
- else if(t<0.0){
- s = clamp(-d/a,0.0,1.0);
- t = 0.0;
- }
- else{
- double invDet = 1.0/det;
- s *= invDet;
- t *= invDet;
- }
+ if(t<0.0){
+ if(d<0.0){
+ s = clamp(-d/a,0.0,1.0);
+ t = 0.0;
+ }
+ else{
+ s = 0.0;
+ t = clamp(-e/c,0.0,1.0);
+ }
+ }
+ else{
+ s = 0.0;
+ t = clamp(-e/c,0.0,1.0);
+ }
+ }
+ else if(t<0.0){
+ s = clamp(-d/a,0.0,1.0);
+ t = 0.0;
+ }
+ else{
+ double invDet = 1.0/det;
+ s *= invDet;
+ t *= invDet;
+ }
}
else{
if(s<0.0){
- double tmp0 = b+d;
- double tmp1 = c+e;
- if(tmp1>tmp0){
- double numer = tmp1-tmp0;
- double denom = a-2.0*b+c;
- s = clamp(numer/denom,0.0,1.0);
- t = 1.0-s;
- }
- else{
- t = clamp(-e/c,0.0,1.0);
- s = 0.0;
- }
- }
- else if(t<0.0){
- if(a+d>b+e){
- double numer = c+e-b-d;
- double denom = a-2.0*b+c;
- s = clamp(numer/denom,0.0,1.0);
- t = 1.0-s;
- }
- else{
- s = clamp(-e/c,0.0,1.0);
- t = 0.0;
- }
- }
- else{
- double numer = c+e-b-d;
- double denom = a-2.0*b+c;
- s = clamp(numer/denom,0.0,1.0);
- t = 1.0-s;
- }
+ double tmp0 = b+d;
+ double tmp1 = c+e;
+ if(tmp1>tmp0){
+ double numer = tmp1-tmp0;
+ double denom = a-2.0*b+c;
+ s = clamp(numer/denom,0.0,1.0);
+ t = 1.0-s;
+ }
+ else{
+ t = clamp(-e/c,0.0,1.0);
+ s = 0.0;
+ }
+ }
+ else if(t<0.0){
+ if(a+d>b+e){
+ double numer = c+e-b-d;
+ double denom = a-2.0*b+c;
+ s = clamp(numer/denom,0.0,1.0);
+ t = 1.0-s;
+ }
+ else{
+ s = clamp(-e/c,0.0,1.0);
+ t = 0.0;
+ }
+ }
+ else{
+ double numer = c+e-b-d;
+ double denom = a-2.0*b+c;
+ s = clamp(numer/denom,0.0,1.0);
+ 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());
+ element->getVertex(0)->z()+s*edge0.z()+t*edge1.z());
}
double Nearest_point::clamp(double x,double min,double max){
@@ -1681,17 +1681,17 @@ void Nearest_point::print_field(GRegion* gr){
for(i=0;i<gr->getNumMeshElements();i++){
element = gr->getMeshElement(i);
- for(j=0;j<element->getNumVertices();j++){
- vertex = element->getVertex(j);
- x = vertex->x();
- y = vertex->y();
- z = vertex->z();
- val = search(x,y,z,vec);
- if(val){
- print_segment(SPoint3(x+k*vec.x(),y+k*vec.y(),z+k*vec.z()),
- SPoint3(x-k*vec.x(),y-k*vec.y(),z-k*vec.z()),file);
- }
- }
+ for(j=0;j<element->getNumVertices();j++){
+ vertex = element->getVertex(j);
+ x = vertex->x();
+ y = vertex->y();
+ z = vertex->z();
+ val = search(x,y,z,vec);
+ if(val){
+ print_segment(SPoint3(x+k*vec.x(),y+k*vec.y(),z+k*vec.z()),
+ SPoint3(x-k*vec.x(),y-k*vec.y(),z-k*vec.z()),file);
+ }
+ }
}
file << "};\n";
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index d5653566b6009e62ac7d32dcc12115fc7d7bf8ae..8b2f3cbd608038dac21841adb9d053df1239c7b1 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -87,29 +87,29 @@ public:
std::list<GEdge*> edges = gf->edges();
std::list<GEdge*>::iterator ite = edges.begin();
while(ite != edges.end()){
- if(!(*ite)->isMeshDegenerated()){
- std::vector<MLine*> temp;
- (*ite)->mesh_vertices.clear();
- for(unsigned int i = 0; i< (*ite)->lines.size(); i+=2){
- if (i+1 >= (*ite)->lines.size()){
- Msg::Error("1D mesh cannot be divided by 2");
- break;
- }
- MVertex *v1 = (*ite)->lines[i]->getVertex(0);
- MVertex *v2 = (*ite)->lines[i]->getVertex(1);
- MVertex *v3 = (*ite)->lines[i+1]->getVertex(1);
- v2->x() = 0.5*(v1->x()+v3->x());
- v2->y() = 0.5*(v1->y()+v3->y());
- v2->z() = 0.5*(v1->z()+v3->z());
- temp.push_back(new MLine(v1,v3));
- if (v1->onWhat() == *ite) (*ite)->mesh_vertices.push_back(v1);
- _middle[MEdge(v1,v3)] = v2;
- }
- _backup[*ite] = (*ite)->lines;
- // printf("line %d goes from %d to %d\n",(*ite)->tag(), (*ite)->lines.size()-1,temp.size()-1);
- (*ite)->lines = temp;
- }
- ++ite;
+ if(!(*ite)->isMeshDegenerated()){
+ std::vector<MLine*> temp;
+ (*ite)->mesh_vertices.clear();
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i+=2){
+ if (i+1 >= (*ite)->lines.size()){
+ Msg::Error("1D mesh cannot be divided by 2");
+ break;
+ }
+ MVertex *v1 = (*ite)->lines[i]->getVertex(0);
+ MVertex *v2 = (*ite)->lines[i]->getVertex(1);
+ MVertex *v3 = (*ite)->lines[i+1]->getVertex(1);
+ v2->x() = 0.5*(v1->x()+v3->x());
+ v2->y() = 0.5*(v1->y()+v3->y());
+ v2->z() = 0.5*(v1->z()+v3->z());
+ temp.push_back(new MLine(v1,v3));
+ if (v1->onWhat() == *ite) (*ite)->mesh_vertices.push_back(v1);
+ _middle[MEdge(v1,v3)] = v2;
+ }
+ _backup[*ite] = (*ite)->lines;
+ // printf("line %d goes from %d to %d\n",(*ite)->tag(), (*ite)->lines.size()-1,temp.size()-1);
+ (*ite)->lines = temp;
+ }
+ ++ite;
}
CTX::instance()->mesh.lcFactor *=2.0;
}
@@ -124,28 +124,28 @@ public:
MVertex *v[3];
SPoint2 m[3];
for (int j=0;j<3;j++){
- MEdge E = _gf->triangles[i]->getEdge(j);
- SPoint2 p1, p2;
- reparamMeshEdgeOnFace(E.getVertex(0),E.getVertex(1),_gf,p1,p2);
- std::map<MEdge, MVertex *, Less_Edge>::iterator it = _middle.find(E);
- std::map<MEdge, MVertex *, Less_Edge>::iterator it2 = eds.find(E);
- m[j] = p1;
- if (it == _middle.end() && it2 == eds.end()){
- GPoint gp = _gf->point((p1+p2)*0.5);
- double XX = 0.5*(E.getVertex(0)->x()+E.getVertex(1)->x());
- double YY = 0.5*(E.getVertex(0)->y()+E.getVertex(1)->y());
- double ZZ = 0.5*(E.getVertex(0)->z()+E.getVertex(1)->z());
- v[j] = new MFaceVertex (XX,YY,ZZ,_gf,gp.u(),gp.v());
- _gf->mesh_vertices.push_back(v[j]);
- eds[E] = v[j];
- }
- else if (it == _middle.end()){
- v[j] = it2->second;
- }
- else {
- v[j] = it->second;
- v[j]->onWhat()->mesh_vertices.push_back(v[j]);
- }
+ MEdge E = _gf->triangles[i]->getEdge(j);
+ SPoint2 p1, p2;
+ reparamMeshEdgeOnFace(E.getVertex(0),E.getVertex(1),_gf,p1,p2);
+ std::map<MEdge, MVertex *, Less_Edge>::iterator it = _middle.find(E);
+ std::map<MEdge, MVertex *, Less_Edge>::iterator it2 = eds.find(E);
+ m[j] = p1;
+ if (it == _middle.end() && it2 == eds.end()){
+ GPoint gp = _gf->point((p1+p2)*0.5);
+ double XX = 0.5*(E.getVertex(0)->x()+E.getVertex(1)->x());
+ double YY = 0.5*(E.getVertex(0)->y()+E.getVertex(1)->y());
+ double ZZ = 0.5*(E.getVertex(0)->z()+E.getVertex(1)->z());
+ v[j] = new MFaceVertex (XX,YY,ZZ,_gf,gp.u(),gp.v());
+ _gf->mesh_vertices.push_back(v[j]);
+ eds[E] = v[j];
+ }
+ else if (it == _middle.end()){
+ v[j] = it2->second;
+ }
+ else {
+ v[j] = it->second;
+ v[j]->onWhat()->mesh_vertices.push_back(v[j]);
+ }
}
GPoint gp = _gf->point((m[0]+m[1]+m[2])*(1./3.));
double XX = (v[0]->x()+v[1]->x()+v[2]->x())*(1./3.);
@@ -163,28 +163,28 @@ public:
MVertex *v[4];
SPoint2 m[4];
for (int j=0;j<4;j++){
- MEdge E = _gf->quadrangles[i]->getEdge(j);
- SPoint2 p1, p2;
- reparamMeshEdgeOnFace(E.getVertex(0),E.getVertex(1),_gf,p1,p2);
- std::map<MEdge, MVertex *, Less_Edge>::iterator it = _middle.find(E);
- std::map<MEdge, MVertex *, Less_Edge>::iterator it2 = eds.find(E);
- m[j] = p1;
- if (it == _middle.end() && it2 == eds.end()){
- GPoint gp = _gf->point((p1+p2)*0.5);
- double XX = 0.5*(E.getVertex(0)->x()+E.getVertex(1)->x());
- double YY = 0.5*(E.getVertex(0)->y()+E.getVertex(1)->y());
- double ZZ = 0.5*(E.getVertex(0)->z()+E.getVertex(1)->z());
- v[j] = new MFaceVertex (XX,YY,ZZ,_gf,gp.u(),gp.v());
- _gf->mesh_vertices.push_back(v[j]);
- eds[E] = v[j];
- }
- else if (it == _middle.end()){
- v[j] = it2->second;
- }
- else {
- v[j] = it->second;
- v[j]->onWhat()->mesh_vertices.push_back(v[j]);
- }
+ MEdge E = _gf->quadrangles[i]->getEdge(j);
+ SPoint2 p1, p2;
+ reparamMeshEdgeOnFace(E.getVertex(0),E.getVertex(1),_gf,p1,p2);
+ std::map<MEdge, MVertex *, Less_Edge>::iterator it = _middle.find(E);
+ std::map<MEdge, MVertex *, Less_Edge>::iterator it2 = eds.find(E);
+ m[j] = p1;
+ if (it == _middle.end() && it2 == eds.end()){
+ GPoint gp = _gf->point((p1+p2)*0.5);
+ double XX = 0.5*(E.getVertex(0)->x()+E.getVertex(1)->x());
+ double YY = 0.5*(E.getVertex(0)->y()+E.getVertex(1)->y());
+ double ZZ = 0.5*(E.getVertex(0)->z()+E.getVertex(1)->z());
+ v[j] = new MFaceVertex (XX,YY,ZZ,_gf,gp.u(),gp.v());
+ _gf->mesh_vertices.push_back(v[j]);
+ eds[E] = v[j];
+ }
+ else if (it == _middle.end()){
+ v[j] = it2->second;
+ }
+ else {
+ v[j] = it->second;
+ v[j]->onWhat()->mesh_vertices.push_back(v[j]);
+ }
}
GPoint gp = _gf->point((m[0]+m[1]+m[2]+m[3])*0.25);
// FIXME : NOT EXACTLY CORRECT, BUT THAT'S THE PLACE WE WANT THE POINT TO RESIDE
@@ -215,11 +215,11 @@ public:
restore();
recombineIntoQuads(_gf,true,true,1.e-3,false);
computeElementShapes(_gf,
- _gf->meshStatistics.worst_element_shape,
- _gf->meshStatistics.average_element_shape,
- _gf->meshStatistics.best_element_shape,
- _gf->meshStatistics.nbTriangle,
- _gf->meshStatistics.nbGoodQuality);
+ _gf->meshStatistics.worst_element_shape,
+ _gf->meshStatistics.average_element_shape,
+ _gf->meshStatistics.best_element_shape,
+ _gf->meshStatistics.nbTriangle,
+ _gf->meshStatistics.nbGoodQuality);
}
}
void restore ()
@@ -228,7 +228,7 @@ public:
std::list<GEdge*>::iterator ite = edges.begin();
while(ite != edges.end()){
for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
- delete (*ite)->lines[i];
+ delete (*ite)->lines[i];
}
(*ite)->lines = _backup[*ite];
++ite;
@@ -354,9 +354,9 @@ static void copyMesh(GFace *source, GFace *target)
}
if (!vt[0] || !vt[1] ||!vt[2]){
Msg::Error("Problem in mesh copying procedure %p %p %p %d %d %d",
- vt[0], vt[1], vt[2], source->triangles[i]->getVertex(0)->onWhat()->dim(),
- source->triangles[i]->getVertex(1)->onWhat()->dim(),
- source->triangles[i]->getVertex(2)->onWhat()->dim());
+ vt[0], vt[1], vt[2], source->triangles[i]->getVertex(0)->onWhat()->dim(),
+ source->triangles[i]->getVertex(1)->onWhat()->dim(),
+ source->triangles[i]->getVertex(2)->onWhat()->dim());
return;
}
target->triangles.push_back(new MTriangle(vt[0], vt[1], vt[2]));
@@ -369,11 +369,11 @@ static void copyMesh(GFace *source, GFace *target)
MVertex *v4 = vs2vt[source->quadrangles[i]->getVertex(3)];
if (!v1 || !v2 || !v3 || !v4){
Msg::Error("Problem in mesh copying procedure %p %p %p %p %d %d %d %d",
- v1, v2, v3, v4,
- source->quadrangles[i]->getVertex(0)->onWhat()->dim(),
- source->quadrangles[i]->getVertex(1)->onWhat()->dim(),
- source->quadrangles[i]->getVertex(2)->onWhat()->dim(),
- source->quadrangles[i]->getVertex(3)->onWhat()->dim());
+ v1, v2, v3, v4,
+ source->quadrangles[i]->getVertex(0)->onWhat()->dim(),
+ source->quadrangles[i]->getVertex(1)->onWhat()->dim(),
+ source->quadrangles[i]->getVertex(2)->onWhat()->dim(),
+ source->quadrangles[i]->getVertex(3)->onWhat()->dim());
}
target->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
}
@@ -549,12 +549,12 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, _fatallyFailed, e2r, notRecovered);
if(e) e->g = g;
else {
- if (_fatallyFailed)
+ if (_fatallyFailed)
Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)",
vstart->x(), vstart->y(), vend->x(), vend->y(), i,
ge->mesh_vertices.size());
- return !_fatallyFailed;
- }
+ return !_fatallyFailed;
+ }
}
}
}
@@ -657,40 +657,40 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
MEdge dv(v1,v2);
addOrRemove(v1,v2,bedges);
for (unsigned int SIDE = 0 ; SIDE < _columns->_normals.count(dv); SIDE ++){
- std::vector<MElement*> myCol;
- edgeColumn ec = _columns->getColumns(v1, v2, SIDE);
- const BoundaryLayerData & c1 = ec._c1;
- const BoundaryLayerData & c2 = ec._c2;
- int N = std::min(c1._column.size(),c2._column.size());
- for (int l=0;l < N ;++l){
- MVertex *v11,*v12,*v21,*v22;
- v21 = c1._column[l];
- v22 = c2._column[l];
- if (l == 0){
- v11 = v1;
- v12 = v2;
- }
- else {
- v11 = c1._column[l-1];
- v12 = c2._column[l-1];
- }
- MEdge dv2 (v21,v22);
-
- //avoid convergent errors
- if (dv2.length() < 0.03 * dv.length())break;
- MQuadrangle *qq = new MQuadrangle(v11,v21,v22,v12);
- myCol.push_back(qq);
- blQuads.push_back(qq);
+ std::vector<MElement*> myCol;
+ edgeColumn ec = _columns->getColumns(v1, v2, SIDE);
+ const BoundaryLayerData & c1 = ec._c1;
+ const BoundaryLayerData & c2 = ec._c2;
+ int N = std::min(c1._column.size(),c2._column.size());
+ for (int l=0;l < N ;++l){
+ MVertex *v11,*v12,*v21,*v22;
+ v21 = c1._column[l];
+ v22 = c2._column[l];
+ if (l == 0){
+ v11 = v1;
+ v12 = v2;
+ }
+ else {
+ v11 = c1._column[l-1];
+ v12 = c2._column[l-1];
+ }
+ MEdge dv2 (v21,v22);
+
+ //avoid convergent errors
+ if (dv2.length() < 0.03 * dv.length())break;
+ MQuadrangle *qq = new MQuadrangle(v11,v21,v22,v12);
+ myCol.push_back(qq);
+ blQuads.push_back(qq);
if(ff2)
fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
v11->x(),v11->y(),v11->z(),
v12->x(),v12->y(),v12->z(),
v22->x(),v22->y(),v22->z(),
v21->x(),v21->y(),v21->z());
- }
- // int M = std::max(c1._column.size(),c2._column.size());
- for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
- _columns->_elemColumns[myCol[0]] = myCol;
+ }
+ // int M = std::max(c1._column.size(),c2._column.size());
+ for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
+ _columns->_elemColumns[myCol[0]] = myCol;
}
}
++ite;
@@ -707,38 +707,38 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
const BoundaryLayerData & c2 = _columns->getColumn(v,i+1);
int N = std::min(c1._column.size(),c2._column.size());
for (int l=0;l < N ;++l){
- MVertex *v11,*v12,*v21,*v22;
- v21 = c1._column[l];
- v22 = c2._column[l];
- if (l == 0){
- v11 = v;
- v12 = v;
- }
- else {
- v11 = c1._column[l-1];
- v12 = c2._column[l-1];
- }
- if (v11 != v12){
- MQuadrangle *qq = new MQuadrangle(v11,v12,v22,v21);
- myCol.push_back(qq);
- blQuads.push_back(qq);
- if(ff2)
+ MVertex *v11,*v12,*v21,*v22;
+ v21 = c1._column[l];
+ v22 = c2._column[l];
+ if (l == 0){
+ v11 = v;
+ v12 = v;
+ }
+ else {
+ v11 = c1._column[l-1];
+ v12 = c2._column[l-1];
+ }
+ if (v11 != v12){
+ MQuadrangle *qq = new MQuadrangle(v11,v12,v22,v21);
+ myCol.push_back(qq);
+ blQuads.push_back(qq);
+ if(ff2)
fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
v11->x(),v11->y(),v11->z(),
v12->x(),v12->y(),v12->z(),
v22->x(),v22->y(),v22->z(),
v21->x(),v21->y(),v21->z());
- }
- else {
- MTriangle *qq = new MTriangle(v,v22,v21);
- myCol.push_back(qq);
- blTris.push_back(qq);
+ }
+ else {
+ MTriangle *qq = new MTriangle(v,v22,v21);
+ myCol.push_back(qq);
+ blTris.push_back(qq);
if(ff2)
fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
v->x(),v->y(),v->z(),
v22->x(),v22->y(),v22->z(),
v21->x(),v21->y(),v21->z());
- }
+ }
}
}
for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
@@ -769,7 +769,7 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
}
discreteEdge ne (gf->model(), 444444,0,
- (*edges.begin())->getEndVertex());
+ (*edges.begin())->getEndVertex());
std::list<GEdge*> hop;
std::set<MEdge,Less_Edge>::iterator it = bedges.begin();
@@ -955,10 +955,10 @@ static void directions_storage(GFace* gf)
// Builds An initial triangular mesh that respects the boundaries of
// the domain, including embedded points and surfaces
bool meshGenerator(GFace *gf, int RECUR_ITER,
- bool repairSelfIntersecting1dMesh,
- bool onlyInitialMesh,
- bool debug,
- std::list<GEdge*> *replacement_edges)
+ bool repairSelfIntersecting1dMesh,
+ bool onlyInitialMesh,
+ bool debug,
+ std::list<GEdge*> *replacement_edges)
{
//onlyInitialMesh=true;
BDS_GeomEntity CLASS_F(1, 2);
@@ -984,11 +984,11 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
if((*ite)->isSeam(gf)) return false;
if(!(*ite)->isMeshDegenerated()){
for(unsigned int i = 0; i< (*ite)->lines.size(); i++)
- medgesToRecover.push_back(MEdge((*ite)->lines[i]->getVertex(0),
- (*ite)->lines[i]->getVertex(1)));
+ medgesToRecover.push_back(MEdge((*ite)->lines[i]->getVertex(0),
+ (*ite)->lines[i]->getVertex(1)));
for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
- MVertex *v1 = (*ite)->lines[i]->getVertex(0);
- MVertex *v2 = (*ite)->lines[i]->getVertex(1);
+ MVertex *v1 = (*ite)->lines[i]->getVertex(0);
+ MVertex *v2 = (*ite)->lines[i]->getVertex(1);
all_vertices.insert(v1);
all_vertices.insert(v2);
if(boundary.find(v1) == boundary.end())
@@ -1018,7 +1018,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
while(ite != emb_edges.end()){
for(unsigned int i = 0; i< (*ite)->lines.size(); i++)
medgesToRecover.push_back(MEdge((*ite)->lines[i]->getVertex(0),
- (*ite)->lines[i]->getVertex(1)));
+ (*ite)->lines[i]->getVertex(1)));
if(!(*ite)->isMeshDegenerated()){
all_vertices.insert((*ite)->mesh_vertices.begin(),
(*ite)->mesh_vertices.end() );
@@ -1055,7 +1055,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
MVertex *vv[3];
int i = 0;
for(std::set<MVertex*, MVertexLessThanNum>::iterator it = all_vertices.begin();
- it != all_vertices.end(); it++){
+ it != all_vertices.end(); it++){
vv[i++] = *it;
}
gf->triangles.push_back(new MTriangle(vv[0], vv[1], vv[2]));
@@ -1102,9 +1102,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
DocRecord doc(points.size() + 4);
for(unsigned int i = 0; i < points.size(); i++){
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
// printf("%22.15E %22.15E \n",XX,YY);
doc.points[i].where.h = points[i]->u + XX;
doc.points[i].where.v = points[i]->v + YY;
@@ -1119,9 +1119,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// add 4 points than encloses the domain (use negative number to
// distinguish those fake vertices)
double bb[4][2] = {{bbox.min().x(), bbox.min().y()},
- {bbox.min().x(), bbox.max().y()},
- {bbox.max().x(), bbox.min().y()},
- {bbox.max().x(), bbox.max().y()}};
+ {bbox.min().x(), bbox.max().y()},
+ {bbox.max().x(), bbox.min().y()},
+ {bbox.max().x(), bbox.max().y()}};
for(int ip = 0; ip < 4; ip++){
BDS_Point *pp = m->add_point(-ip - 1, bb[ip][0], bb[ip][1], gf);
m->add_geom(gf->tag(), 2);
@@ -1156,8 +1156,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
int c = doc.triangles[i].c;
int n = doc.numPoints;
if(a < 0 || a >= n || b < 0 || b >= n || c < 0 || c >= n){
- Msg::Warning("Skipping bad triangle %d", i);
- continue;
+ Msg::Warning("Skipping bad triangle %d", i);
+ continue;
}
BDS_Point *p1 = (BDS_Point*)doc.points[doc.triangles[i].a].data;
BDS_Point *p2 = (BDS_Point*)doc.points[doc.triangles[i].b].data;
@@ -1242,31 +1242,31 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
while(ite != edges.end()){
if(!(*ite)->isMeshDegenerated()){
if (!recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2)){
- delete m;
- gf->meshStatistics.status = GFace::FAILED;
- return false;
+ delete m;
+ gf->meshStatistics.status = GFace::FAILED;
+ return false;
}
}
++ite;
}
Msg::Debug("Recovering %d mesh Edges (%d not recovered)", edgesToRecover.size(),
- edgesNotRecovered.size());
+ edgesNotRecovered.size());
if(edgesNotRecovered.size()){
std::ostringstream sstream;
for(std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
- itr != edgesNotRecovered.end(); ++itr)
+ itr != edgesNotRecovered.end(); ++itr)
sstream << " " << itr->ge->tag();
Msg::Warning(":-( There are %d intersections in the 1D mesh (curves%s)",
- edgesNotRecovered.size(), sstream.str().c_str());
+ edgesNotRecovered.size(), sstream.str().c_str());
if (repairSelfIntersecting1dMesh)
Msg::Warning("8-| Gmsh splits those edges and tries again");
if(debug){
char name[245];
sprintf(name, "surface%d-not_yet_recovered-real-%d.msh", gf->tag(),
- RECUR_ITER);
+ RECUR_ITER);
gf->model()->writeMSH(name);
}
@@ -1278,14 +1278,14 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
//int *_error = new int[3 * edgesNotRecovered.size()];
int I = 0;
for(; itr != edgesNotRecovered.end(); ++itr){
- int p1 = itr->p1;
- int p2 = itr->p2;
- int tag = itr->ge->tag();
- Msg::Error("Edge not recovered: %d %d %d", p1, p2, tag);
- //_error[3 * I + 0] = p1;
- //_error[3 * I + 1] = p2;
- //_error[3 * I + 2] = tag;
- I++;
+ int p1 = itr->p1;
+ int p2 = itr->p2;
+ int tag = itr->ge->tag();
+ Msg::Error("Edge not recovered: %d %d %d", p1, p2, tag);
+ //_error[3 * I + 0] = p1;
+ //_error[3 * I + 1] = p2;
+ //_error[3 * I + 2] = tag;
+ I++;
}
//throw _error;
}
@@ -1294,14 +1294,14 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
delete m;
if(RECUR_ITER < 10 && facesToRemesh.size() == 0)
return meshGenerator
- (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh,
- debug, replacement_edges);
+ (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh,
+ debug, replacement_edges);
return false;
}
if(RECUR_ITER > 0)
Msg::Warning(":-) Gmsh was able to recover all edges after %d iterations",
- RECUR_ITER);
+ RECUR_ITER);
Msg::Debug("Boundary Edges recovered for surface %d", gf->tag());
@@ -1319,9 +1319,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
BDS_Point *n[4];
t->getNodes(n);
if (n[0]->iD < 0 || n[1]->iD < 0 ||
- n[2]->iD < 0 ) {
- recur_tag(t, &CLASS_EXTERIOR);
- break;
+ n[2]->iD < 0 ) {
+ recur_tag(t, &CLASS_EXTERIOR);
+ break;
}
++itt;
}
@@ -1334,14 +1334,14 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
while (ite != m->edges.end()){
BDS_Edge *e = *ite;
if(e->g && e->numfaces() == 2){
- if(e->faces(0)->g == &CLASS_EXTERIOR){
- recur_tag(e->faces(1), &CLASS_F);
- break;
- }
- else if(e->faces(1)->g == &CLASS_EXTERIOR){
- recur_tag(e->faces(0), &CLASS_F);
- break;
- }
+ if(e->faces(0)->g == &CLASS_EXTERIOR){
+ recur_tag(e->faces(1), &CLASS_F);
+ break;
+ }
+ else if(e->faces(1)->g == &CLASS_EXTERIOR){
+ recur_tag(e->faces(0), &CLASS_F);
+ break;
+ }
}
++ite;
}
@@ -1357,16 +1357,16 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
while (ite != m->edges.end()){
BDS_Edge *e = *ite;
if(e->g && e->numfaces() == 2){
- BDS_Point *oface[2];
- e->oppositeof(oface);
- if(oface[0]->iD < 0){
- recur_tag(e->faces(1), &CLASS_F);
- break;
- }
- else if(oface[1]->iD < 0){
- recur_tag(e->faces(0), &CLASS_F);
- break;
- }
+ BDS_Point *oface[2];
+ e->oppositeof(oface);
+ if(oface[0]->iD < 0){
+ recur_tag(e->faces(1), &CLASS_F);
+ break;
+ }
+ else if(oface[1]->iD < 0){
+ recur_tag(e->faces(0), &CLASS_F);
+ break;
+ }
}
++ite;
}
@@ -1382,28 +1382,28 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// compute characteristic lengths at vertices
if (!onlyInitialMesh){
Msg::Debug("Computing mesh size field at mesh vertices %d",
- edgesToRecover.size());
+ edgesToRecover.size());
std::set<BDS_Point*, PointLessThan>::iterator it = m->points.begin();
for(; it != m->points.end();++it){
- // for(int i = 0; i < doc.numPoints; i++){
- // BDS_Point *pp = (BDS_Point*)doc.points[i].data;
- BDS_Point *pp = *it;
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itv = recoverMap.find(pp);
- if(itv != recoverMap.end()){
- MVertex *here = itv->second;
- GEntity *ge = here->onWhat();
- if(ge->dim() == 0){
- pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
- }
- else if(ge->dim() == 1){
- double u;
- here->getParameter(0, u);
- pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
- }
- else
- pp->lcBGM() = MAX_LC;
- pp->lc() = pp->lcBGM();
- }
+ // for(int i = 0; i < doc.numPoints; i++){
+ // BDS_Point *pp = (BDS_Point*)doc.points[i].data;
+ BDS_Point *pp = *it;
+ std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itv = recoverMap.find(pp);
+ if(itv != recoverMap.end()){
+ MVertex *here = itv->second;
+ GEntity *ge = here->onWhat();
+ if(ge->dim() == 0){
+ pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
+ }
+ else if(ge->dim() == 1){
+ double u;
+ here->getParameter(0, u);
+ pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
+ }
+ else
+ pp->lcBGM() = MAX_LC;
+ pp->lc() = pp->lcBGM();
+ }
}
}
@@ -1503,10 +1503,10 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
/*
computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index,
- gf->meshStatistics.longest_edge_length,
- gf->meshStatistics.smallest_edge_length,
- gf->meshStatistics.nbEdge,
- gf->meshStatistics.nbGoodLength);
+ gf->meshStatistics.longest_edge_length,
+ gf->meshStatistics.smallest_edge_length,
+ gf->meshStatistics.nbEdge,
+ gf->meshStatistics.nbGoodLength);
*/
//printf("=== Efficiency index is tau=%g\n", gf->meshStatistics.efficiency_index);
@@ -1575,8 +1575,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
gf->meshStatistics.nbGoodQuality);
gf->mesh_vertices.insert(gf->mesh_vertices.end(),
- gf->additionalVertices.begin(),
- gf->additionalVertices.end());
+ gf->additionalVertices.begin(),
+ gf->additionalVertices.end());
gf->additionalVertices.clear();
if(CTX::instance()->mesh.algo3d==ALGO_3D_RTREE){
@@ -1925,7 +1925,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
MVertex *vv[3];
int i = 0;
for(std::map<BDS_Point*, MVertex*, PointLessThan>::iterator it = recoverMap.begin();
- it != recoverMap.end(); it++){
+ it != recoverMap.end(); it++){
vv[i++] = it->second;
}
gf->triangles.push_back(new MTriangle(vv[0], vv[1], vv[2]));
@@ -2022,8 +2022,8 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
BDS_Point *pp = edgeLoop_BDS[j];
- v.push_back(recoverMap[pp]);
- recoverMapInv[recoverMap[pp]] = pp;
+ v.push_back(recoverMap[pp]);
+ recoverMapInv[recoverMap[pp]] = pp;
}
}
@@ -2237,29 +2237,29 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
BDS_Point *bds = it->first;
std::map<MVertex*, BDS_Point*>::iterator invIt = invertMap.find(mv1);
if (invIt != invertMap.end()){
- // create a new "fake" vertex that will be destroyed afterwards
- MVertex *mv2 = 0;
- if (mv1->onWhat()->dim() == 1) {
- double t;
- mv1->getParameter(0,t);
- mv2 = new MEdgeVertex(mv1->x(),mv1->y(),mv1->z(),mv1->onWhat(), t,
- ((MEdgeVertex*)mv1)->getLc());
- }
- else if (mv1->onWhat()->dim() == 0) {
- mv2 = new MVertex (mv1->x(),mv1->y(),mv1->z(),mv1->onWhat());
- }
- else
- Msg::Error("Could not reconstruct seam");
- if(mv2){
- it->second = mv2;
- equivalence[mv2] = mv1;
- parametricCoordinates[mv2] = SPoint2(bds->u,bds->v);
- invertMap[mv2] = bds;
- }
+ // create a new "fake" vertex that will be destroyed afterwards
+ MVertex *mv2 = 0;
+ if (mv1->onWhat()->dim() == 1) {
+ double t;
+ mv1->getParameter(0,t);
+ mv2 = new MEdgeVertex(mv1->x(),mv1->y(),mv1->z(),mv1->onWhat(), t,
+ ((MEdgeVertex*)mv1)->getLc());
+ }
+ else if (mv1->onWhat()->dim() == 0) {
+ mv2 = new MVertex (mv1->x(),mv1->y(),mv1->z(),mv1->onWhat());
+ }
+ else
+ Msg::Error("Could not reconstruct seam");
+ if(mv2){
+ it->second = mv2;
+ equivalence[mv2] = mv1;
+ parametricCoordinates[mv2] = SPoint2(bds->u,bds->v);
+ invertMap[mv2] = bds;
+ }
}
else {
- parametricCoordinates[mv1] = SPoint2(bds->u,bds->v);
- invertMap[mv1] = bds;
+ parametricCoordinates[mv1] = SPoint2(bds->u,bds->v);
+ invertMap[mv1] = bds;
}
++it;
}
@@ -2298,11 +2298,11 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// may be created, for example on a sphere that contains one
// pole
if(v1 != v2 && v1 != v3 && v2 != v3){
- // we are in the periodic case. if we aim at
- // using delaunay mesh generation in thoses cases,
- // we should double some of the vertices
+ // we are in the periodic case. if we aim at
+ // using delaunay mesh generation in thoses cases,
+ // we should double some of the vertices
gf->triangles.push_back(new MTriangle(v1, v2, v3));
- }
+ }
}
else{
MVertex *v4 = recoverMap[n[3]];
@@ -2396,6 +2396,8 @@ void meshGFace::operator() (GFace *gf, bool print)
deMeshGFace dem;
dem(gf);
+ // FIXME: if transfinite surface, impossible to use ALGO_3D_RTREE
+ // because meshGenerator never called
if(MeshTransfiniteSurface(gf)) return;
if(MeshExtrudedSurface(gf)) return;
if(gf->meshMaster() != gf){
@@ -2455,7 +2457,7 @@ void meshGFace::operator() (GFace *gf, bool print)
(noSeam(gf) || gf->getNativeType() == GEntity::GmshModel ||
gf->edgeLoops.empty())){
meshGenerator(gf, 0, repairSelfIntersecting1dMesh, onlyInitialMesh,
- debugSurface >= 0 || debugSurface == -100);
+ debugSurface >= 0 || debugSurface == -100);
}
else {
if(!meshGeneratorPeriodic
@@ -2524,7 +2526,7 @@ void partitionAndRemesh(GFaceCompound *gf)
int allowType = gf->allowPartition();
multiscalePartition *msp = new multiscalePartition(elements, abs(gf->nbSplit),
- method, allowType);
+ method, allowType);
int NF = msp->getNumberOfParts();
int numv = gf->model()->getMaxElementaryNumber(0) + 1;
@@ -2569,7 +2571,7 @@ void partitionAndRemesh(GFaceCompound *gf)
num_gfc, pf->tag());
GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound, U0,
- gf->getTypeOfCompound());
+ gf->getTypeOfCompound());
gfc->meshAttributes.recombine = gf->meshAttributes.recombine;
gf->model()->add(gfc);
@@ -2617,9 +2619,9 @@ void partitionAndRemesh(GFaceCompound *gf)
//update mesh statistics
gf->meshStatistics.efficiency_index += gfc->meshStatistics.efficiency_index;
gf->meshStatistics.longest_edge_length = std::max(gf->meshStatistics.longest_edge_length,
- gfc->meshStatistics.longest_edge_length);
+ gfc->meshStatistics.longest_edge_length);
gf->meshStatistics.smallest_edge_length= std::min(gf->meshStatistics.smallest_edge_length,
- gfc->meshStatistics.smallest_edge_length);
+ gfc->meshStatistics.smallest_edge_length);
gf->meshStatistics.nbGoodLength += gfc->meshStatistics.nbGoodLength;
gf->meshStatistics.nbGoodQuality += gfc->meshStatistics.nbGoodQuality;
gf->meshStatistics.nbEdge += gfc->meshStatistics.nbEdge;
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index ed134404a4ddddcfed278ed791fad1529d7bbdb0..ef22671c94205569b4954f811ef1c15358c2dc54 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -92,14 +92,14 @@ static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes, bidimMeshDa
MVertex *vi = t->getVertex(i);
MVertex *vj = t->getVertex(j);
if (vi != data.equivalent(vj) && vj != data.equivalent(vi) ){
- double dx = vi->x() - vj->x();
- double dy = vi->y() - vj->y();
- double dz = vi->z() - vj->z();
- double l = sqrt(dx * dx + dy * dy + dz * dz);
- std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
- std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
- if(iti->second < 0 || iti->second > l) iti->second = l;
- if(itj->second < 0 || itj->second > l) itj->second = l;
+ double dx = vi->x() - vj->x();
+ double dy = vi->y() - vj->y();
+ double dz = vi->z() - vj->z();
+ double l = sqrt(dx * dx + dy * dy + dz * dz);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ if(iti->second < 0 || iti->second > l) iti->second = l;
+ if(itj->second < 0 || itj->second > l) itj->second = l;
}
}
}
@@ -107,7 +107,7 @@ static void setLcs(MTriangle *t, std::map<MVertex*, double> &vSizes, bidimMeshDa
void buildMeshGenerationDataStructures(GFace *gf,
std::set<MTri3*, compareTri3Ptr> &AllTris,
- bidimMeshData & data)
+ bidimMeshData & data)
{
std::map<MVertex*, double> vSizesMap;
std::list<GEdge*> edges = gf->edges();
@@ -140,9 +140,9 @@ void buildMeshGenerationDataStructures(GFace *gf,
std::list<GEdge*>::iterator ite = embedded_edges.begin();
while(ite != embedded_edges.end()){
if(!(*ite)->isMeshDegenerated()){
- for (unsigned int i = 0; i < (*ite)->lines.size(); i++)
- data.internalEdges.insert(MEdge((*ite)->lines[i]->getVertex(0),
- (*ite)->lines[i]->getVertex(1)));
+ for (unsigned int i = 0; i < (*ite)->lines.size(); i++)
+ data.internalEdges.insert(MEdge((*ite)->lines[i]->getVertex(0),
+ (*ite)->lines[i]->getVertex(1)));
}
++ite;
}
@@ -173,11 +173,11 @@ void computeEquivalences(GFace *gf, bidimMeshData & data)
MTriangle *t = gf->triangles[i];
MVertex *v[3];
for (int j=0;j<3;j++){
- v[j] = t->getVertex(j);
- std::map<MVertex* , MVertex*>::iterator it = data.equivalence->find(v[j]);
- if (it != data.equivalence->end()){
- v[j] = it->second;
- }
+ v[j] = t->getVertex(j);
+ std::map<MVertex* , MVertex*>::iterator it = data.equivalence->find(v[j]);
+ if (it != data.equivalence->end()){
+ v[j] = it->second;
+ }
}
if (v[0] != v[1] && v[0] != v[2] && v[2] != v[1])
newT.push_back(new MTriangle (v[0],v[1],v[2]));
@@ -210,7 +210,7 @@ struct equivalentTriangle {
};
bool computeEquivalentTriangles (GFace *gf,
- std::map<MVertex* , MVertex*>* equivalence)
+ std::map<MVertex* , MVertex*>* equivalence)
{
if (!equivalence)return false;
std::vector<MTriangle*> WTF;
@@ -277,8 +277,8 @@ void transferDataStructure(GFace *gf, std::set<MTri3*, compareTri3Ptr> &AllTris,
index1 = data.getIndex (v1);
index2 = data.getIndex (v2);
normal3points(data.Us[index0], data.Vs[index0], 0.,
- data.Us[index1], data.Vs[index1], 0.,
- data.Us[index2], data.Vs[index2], 0., n2);
+ data.Us[index1], data.Vs[index1], 0.,
+ data.Us[index2], data.Vs[index2], 0., n2);
double pp; prosca(n1, n2, &pp);
if(pp < 0) t->reverse();
}
@@ -420,7 +420,7 @@ int _removeThreeTrianglesNodes(GFace *gf)
n++;
gf->triangles.push_back(newt);
for(int i=0;i<3;i++) {
- touched.insert(lt[i]);
+ touched.insert(lt[i]);
}
}
it++;
@@ -484,18 +484,18 @@ static int _removeTwoQuadsNodes(GFace *gf)
q2->writePOS(stdout,true,false,false,false,false,false);
return 0;
}
- MQuadrangle *q = new MQuadrangle(v1,v2,v3,v4);
- double s1 = surfaceFaceUV(q,gf);
- double s2 = surfaceFaceUV(q1,gf) + surfaceFaceUV(q2,gf);;
- if (s1 > s2){
- delete q;
- }
- else{
- touched.insert(q1);
- touched.insert(q2);
- gf->quadrangles.push_back(q);
- vtouched.insert(it->first);
- }
+ MQuadrangle *q = new MQuadrangle(v1,v2,v3,v4);
+ double s1 = surfaceFaceUV(q,gf);
+ double s2 = surfaceFaceUV(q1,gf) + surfaceFaceUV(q2,gf);;
+ if (s1 > s2){
+ delete q;
+ }
+ else{
+ touched.insert(q1);
+ touched.insert(q2);
+ gf->quadrangles.push_back(q);
+ vtouched.insert(it->first);
+ }
}
}
it++;
@@ -539,11 +539,11 @@ int removeTwoQuadsNodes(GFace *gf)
// collapse v1 & v2 to their middle and replace into e1 & e2
static bool _tryToCollapseThatVertex (GFace *gf,
- std::vector<MElement*> &e1,
- std::vector<MElement*> &e2,
- MElement *q,
- MVertex *v1,
- MVertex *v2)
+ std::vector<MElement*> &e1,
+ std::vector<MElement*> &e2,
+ MElement *q,
+ MVertex *v1,
+ MVertex *v2)
{
std::vector<MElement*> e = e1;
@@ -599,11 +599,11 @@ static bool _tryToCollapseThatVertex (GFace *gf,
v1->setParameter(0,pp.u());v1->setParameter(1,pp.v());
for (unsigned int j=0;j<e.size();++j){
if (e[j] != q){
- for (int k=0;k<4;k++){
- if (e[j]->getVertex(k) == v2){
- e[j]->setVertex(k,v1);
- }
- }
+ for (int k=0;k<4;k++){
+ if (e[j]->getVertex(k) == v2){
+ e[j]->setVertex(k,v1);
+ }
+ }
}
}
return true;
@@ -683,32 +683,32 @@ static int _removeDiamonds(GFace *gf)
touched.find(v3) == touched.end() &&
touched.find(v4) == touched.end() ) {
if (v1->onWhat()->dim() == 2 &&
- v2->onWhat()->dim() == 2 &&
- v3->onWhat()->dim() == 2 &&
- v4->onWhat()->dim() == 2 &&
- it1->second.size() ==3 && it3->second.size() == 3 &&
- _tryToCollapseThatVertex (gf, it1->second, it3->second,
- q, v1, v3)){
- touched.insert(v1);
- touched.insert(v2);
- touched.insert(v3);
- touched.insert(v4);
- deleted.insert(v3);
- diamonds.insert(q);
+ v2->onWhat()->dim() == 2 &&
+ v3->onWhat()->dim() == 2 &&
+ v4->onWhat()->dim() == 2 &&
+ it1->second.size() ==3 && it3->second.size() == 3 &&
+ _tryToCollapseThatVertex (gf, it1->second, it3->second,
+ q, v1, v3)){
+ touched.insert(v1);
+ touched.insert(v2);
+ touched.insert(v3);
+ touched.insert(v4);
+ deleted.insert(v3);
+ diamonds.insert(q);
}
else if (v1->onWhat()->dim() == 2 &&
- v2->onWhat()->dim() == 2 &&
- v3->onWhat()->dim() == 2 &&
- v4->onWhat()->dim() == 2 &&
- it2->second.size() ==3 && it4->second.size() == 3 &&
- _tryToCollapseThatVertex (gf, it2->second, it4->second,
- q, v2, v4)){
- touched.insert(v1);
- touched.insert(v2);
- touched.insert(v3);
- touched.insert(v4);
- deleted.insert(v4);
- diamonds.insert(q);
+ v2->onWhat()->dim() == 2 &&
+ v3->onWhat()->dim() == 2 &&
+ v4->onWhat()->dim() == 2 &&
+ it2->second.size() ==3 && it4->second.size() == 3 &&
+ _tryToCollapseThatVertex (gf, it2->second, it4->second,
+ q, v2, v4)){
+ touched.insert(v1);
+ touched.insert(v2);
+ touched.insert(v3);
+ touched.insert(v4);
+ deleted.insert(v4);
+ diamonds.insert(q);
}
else {
quadrangles2.push_back(q);
@@ -770,12 +770,12 @@ void _relocateVertex(GFace *gf, MVertex *ver,
MEdge e = lt[i]->getEdge(j);
SPoint2 param0, param1;
if (e.getVertex(0) == ver){
- reparamMeshEdgeOnFace(e.getVertex(0), e.getVertex(1), gf, param0, param1);
- pts[e.getVertex(1)] = param1;
+ reparamMeshEdgeOnFace(e.getVertex(0), e.getVertex(1), gf, param0, param1);
+ pts[e.getVertex(1)] = param1;
}
else if (e.getVertex(1) == ver){
- reparamMeshEdgeOnFace(e.getVertex(0), e.getVertex(1), gf, param0, param1);
- pts[e.getVertex(0)] = param0;
+ reparamMeshEdgeOnFace(e.getVertex(0), e.getVertex(1), gf, param0, param1);
+ pts[e.getVertex(0)] = param0;
}
}
}
@@ -791,8 +791,8 @@ void _relocateVertex(GFace *gf, MVertex *ver,
d.normalize();
buildMetric (gf,adj,metric);
const double F = sqrt(metric[0]*d.x()*d.x() +
- 2*metric[1]*d.x()*d.y() +
- metric[2]*d.y()*d.y());
+ 2*metric[1]*d.x()*d.y() +
+ metric[2]*d.y()*d.y());
// printf("%g ",F);
after += adj*F;
COUNT += F;
@@ -809,12 +809,12 @@ void _relocateVertex(GFace *gf, MVertex *ver,
if (success){
GPoint pt = gf->point(trial);
if(pt.succeeded()){
- actual = trial;
- ver->setParameter(0, trial.x());
- ver->setParameter(1, trial.y());
- ver->x() = pt.x();
- ver->y() = pt.y();
- ver->z() = pt.z();
+ actual = trial;
+ ver->setParameter(0, trial.x());
+ ver->setParameter(1, trial.y());
+ ver->x() = pt.x();
+ ver->y() = pt.y();
+ ver->z() = pt.z();
}
}
FACTOR /= 1.4;
@@ -888,11 +888,11 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
const double triQuality = std::min(qmTriangle::gamma(t1b),
qmTriangle::gamma(t2b));
if (!edgeSwapDelProj(v1,v2,v3,v4)){
- if(triQuality < triQualityRef){
- delete t1b;
- delete t2b;
- return false;
- }
+ if(triQuality < triQualityRef){
+ delete t1b;
+ delete t2b;
+ return false;
+ }
}
break;
}
@@ -1120,8 +1120,8 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
if (pairs[i].n2->onWhat()->dim() < 2) NB++;
if (pairs[i].n3->onWhat()->dim() < 2) NB++;
if (pairs[i].n4->onWhat()->dim() < 2) NB++;
- if (elen[i] > (int)1000*exp(.1) && NB > 2) { elen[i] = 5000; }
- else if (elen[i] >= 1000 && NB > 2) { elen[i] = 10000; }
+ if (elen[i] > (int)1000*exp(.1) && NB > 2) { elen[i] = 5000; }
+ else if (elen[i] >= 1000 && NB > 2) { elen[i] = 10000; }
}
if (cubicGraph){
@@ -1140,9 +1140,9 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
sprintf(MATCHFILE,".face.match");
if(perfect_match(ncount, NULL, ecount, &elist, &elen, NULL, MATCHFILE,
0, 0, 0, 0, &matzeit)){
- Msg::Error("Perfect Match failed in Quadrangulation, try something else");
- free(elist);
- pairs.clear();
+ Msg::Error("Perfect Match failed in Quadrangulation, try something else");
+ free(elist);
+ pairs.clear();
}
else{
// TEST
@@ -1150,8 +1150,8 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
int i1 = elist[1+3*k], i2 = elist[1+3*k+1], an=elist[1+3*k+2];
// FIXME !
if (an == 100000 /*|| an == 1000*/){
- // toProcess.push_back(std::make_pair(n2t[i1],n2t[i2]));
- // Msg::Warning("Extra edge found in blossom algorithm, optimization will be required");
+ // toProcess.push_back(std::make_pair(n2t[i1],n2t[i2]));
+ // Msg::Warning("Extra edge found in blossom algorithm, optimization will be required");
}
else{
MElement *t1 = n2t[i1];
@@ -1185,7 +1185,7 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
}
free(elist);
pairs.clear();
- Msg::Debug("Perfect Match Succeeded in Quadrangulation (%g sec)", matzeit);
+ Msg::Debug("Perfect Match Succeeded in Quadrangulation (%g sec)", matzeit);
}
}
@@ -1269,8 +1269,8 @@ static double printStats(GFace *gf,const char *message)
void recombineIntoQuads(GFace *gf,
bool topologicalOpti,
bool nodeRepositioning,
- double minqual,
- bool firstpass)
+ double minqual,
+ bool firstpass)
{
double t1 = Cpu();
@@ -1294,15 +1294,15 @@ void recombineIntoQuads(GFace *gf,
if (saveAll) gf->model()->writeMSH("smoothed.msh");
int ITER=0;
std::set<MEdge,Less_Edge> prioritory;
- int nbTwoQuadNodes = 1;
- int nbDiamonds = 1;
- while(nbTwoQuadNodes || nbDiamonds){
- nbTwoQuadNodes = removeTwoQuadsNodes(gf);
- nbDiamonds = removeDiamonds(gf) ;
- if(haveParam) RelocateVertices (gf,CTX::instance()->mesh.nbSmoothing);
- // printStats (gf, "toto");
- if (ITER > 20) break;
- ITER ++;
+ int nbTwoQuadNodes = 1;
+ int nbDiamonds = 1;
+ while(nbTwoQuadNodes || nbDiamonds){
+ nbTwoQuadNodes = removeTwoQuadsNodes(gf);
+ nbDiamonds = removeDiamonds(gf) ;
+ if(haveParam) RelocateVertices (gf,CTX::instance()->mesh.nbSmoothing);
+ // printStats (gf, "toto");
+ if (ITER > 20) break;
+ ITER ++;
}
}
}
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 2d1bc52f6148dcad6c2c5436eaf6320a5f538baa..106361cc392ce06c4d80aec39620756035a3bcd1 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -70,48 +70,48 @@ class splitQuadRecovery {
for (std::multimap<GEntity*, std::pair<MVertex*,MFace> >::iterator it2 =
_data.lower_bound(*it); it2 != _data.upper_bound(*it) ; ++it2){
const MFace &f = it2->second.second;
- MVertex *v = it2->second.first;
- v->onWhat()->mesh_vertices.erase(std::find(v->onWhat()->mesh_vertices.begin(),
- v->onWhat()->mesh_vertices.end(), v));
- std::set<MFace, Less_Face>::iterator itf0 = allFaces.find(MFace(f.getVertex(0),
- f.getVertex(1),v));
- std::set<MFace, Less_Face>::iterator itf1 = allFaces.find(MFace(f.getVertex(1),
- f.getVertex(2),v));
- std::set<MFace, Less_Face>::iterator itf2 = allFaces.find(MFace(f.getVertex(2),
- f.getVertex(3),v));
- std::set<MFace, Less_Face>::iterator itf3 = allFaces.find(MFace(f.getVertex(3),
- f.getVertex(0),v));
- if (itf0 != allFaces.end() && itf1 != allFaces.end() &&
- itf2 != allFaces.end() && itf3 != allFaces.end()){
- (*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0), f.getVertex(1),
- f.getVertex(2), f.getVertex(3)));
- allFaces.erase(*itf0);
- allFaces.erase(*itf1);
- allFaces.erase(*itf2);
- allFaces.erase(*itf3);
- // printf("some pyramids should be created %d regions\n", (*it)->numRegions());
- for (int iReg = 0; iReg < (*it)->numRegions(); iReg++){
- if (iReg == 1) {
- Msg::Error("Cannot build pyramids on non manifold faces");
- v = new MVertex(v->x(), v->y(), v->z(), (*it)->getRegion(iReg));
- }
- else
- v->setEntity ((*it)->getRegion(iReg));
- // A quad face connected to an hex or a primsm --> leave the quad face as is
- if (_toDelete.find(f) == _toDelete.end()){
- (*it)->getRegion(iReg)->pyramids.push_back
- (new MPyramid(f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), v));
- (*it)->getRegion(iReg)->mesh_vertices.push_back(v);
- NBPY++;
- }
- else {
- delete v;
- }
- }
- }
+ MVertex *v = it2->second.first;
+ v->onWhat()->mesh_vertices.erase(std::find(v->onWhat()->mesh_vertices.begin(),
+ v->onWhat()->mesh_vertices.end(), v));
+ std::set<MFace, Less_Face>::iterator itf0 = allFaces.find(MFace(f.getVertex(0),
+ f.getVertex(1),v));
+ std::set<MFace, Less_Face>::iterator itf1 = allFaces.find(MFace(f.getVertex(1),
+ f.getVertex(2),v));
+ std::set<MFace, Less_Face>::iterator itf2 = allFaces.find(MFace(f.getVertex(2),
+ f.getVertex(3),v));
+ std::set<MFace, Less_Face>::iterator itf3 = allFaces.find(MFace(f.getVertex(3),
+ f.getVertex(0),v));
+ if (itf0 != allFaces.end() && itf1 != allFaces.end() &&
+ itf2 != allFaces.end() && itf3 != allFaces.end()){
+ (*it)->quadrangles.push_back(new MQuadrangle(f.getVertex(0), f.getVertex(1),
+ f.getVertex(2), f.getVertex(3)));
+ allFaces.erase(*itf0);
+ allFaces.erase(*itf1);
+ allFaces.erase(*itf2);
+ allFaces.erase(*itf3);
+ // printf("some pyramids should be created %d regions\n", (*it)->numRegions());
+ for (int iReg = 0; iReg < (*it)->numRegions(); iReg++){
+ if (iReg == 1) {
+ Msg::Error("Cannot build pyramids on non manifold faces");
+ v = new MVertex(v->x(), v->y(), v->z(), (*it)->getRegion(iReg));
+ }
+ else
+ v->setEntity ((*it)->getRegion(iReg));
+ // A quad face connected to an hex or a primsm --> leave the quad face as is
+ if (_toDelete.find(f) == _toDelete.end()){
+ (*it)->getRegion(iReg)->pyramids.push_back
+ (new MPyramid(f.getVertex(0), f.getVertex(1), f.getVertex(2), f.getVertex(3), v));
+ (*it)->getRegion(iReg)->mesh_vertices.push_back(v);
+ NBPY++;
+ }
+ else {
+ delete v;
+ }
+ }
+ }
}
for (std::set<MFace, Less_Face>::iterator itf = allFaces.begin();
- itf != allFaces.end(); ++itf){
+ itf != allFaces.end(); ++itf){
(*it)->triangles.push_back
(new MTriangle(itf->getVertex(0), itf->getVertex(1), itf->getVertex(2)));
}
@@ -334,76 +334,76 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
double guess[2] = {0, 0};
if (needParam) {
- int Count = 0;
- for(int j = 0; j < 3; j++){
- if(!v[j]->onWhat()){
- Msg::Error("Uncategorized vertex %d", v[j]->getNum());
- }
- else if(v[j]->onWhat()->dim() == 2){
- double uu,vv;
- v[j]->getParameter(0, uu);
- v[j]->getParameter(1, vv);
- guess[0] += uu;
- guess[1] += vv;
- Count++;
- }
- else if(v[j]->onWhat()->dim() == 1){
- GEdge *ge = (GEdge*)v[j]->onWhat();
- double UU;
- v[j]->getParameter(0, UU);
- SPoint2 param;
- param = ge->reparamOnFace(gf, UU, 1);
- guess[0] += param.x();
- guess[1] += param.y();
- Count++;
- }
- else if(v[j]->onWhat()->dim() == 0){
- SPoint2 param;
- GVertex *gv = (GVertex*)v[j]->onWhat();
- param = gv->reparamOnFace(gf,1);
- guess[0] += param.x();
- guess[1] += param.y();
- Count++;
- }
- }
- if(Count != 0){
- guess[0] /= Count;
- guess[1] /= Count;
- }
+ int Count = 0;
+ for(int j = 0; j < 3; j++){
+ if(!v[j]->onWhat()){
+ Msg::Error("Uncategorized vertex %d", v[j]->getNum());
+ }
+ else if(v[j]->onWhat()->dim() == 2){
+ double uu,vv;
+ v[j]->getParameter(0, uu);
+ v[j]->getParameter(1, vv);
+ guess[0] += uu;
+ guess[1] += vv;
+ Count++;
+ }
+ else if(v[j]->onWhat()->dim() == 1){
+ GEdge *ge = (GEdge*)v[j]->onWhat();
+ double UU;
+ v[j]->getParameter(0, UU);
+ SPoint2 param;
+ param = ge->reparamOnFace(gf, UU, 1);
+ guess[0] += param.x();
+ guess[1] += param.y();
+ Count++;
+ }
+ else if(v[j]->onWhat()->dim() == 0){
+ SPoint2 param;
+ GVertex *gv = (GVertex*)v[j]->onWhat();
+ param = gv->reparamOnFace(gf,1);
+ guess[0] += param.x();
+ guess[1] += param.y();
+ Count++;
+ }
+ }
+ if(Count != 0){
+ guess[0] /= Count;
+ guess[1] /= Count;
+ }
}
for(int j = 0; j < 3; j++){
- if (out.trifacelist[i * 3 + j] - 1 >= initialSize){
- printf("aaaaaaaaaaaaaaargh\n");
- // if(v[j]->onWhat()->dim() == 3){
- v[j]->onWhat()->mesh_vertices.erase
- (std::find(v[j]->onWhat()->mesh_vertices.begin(),
- v[j]->onWhat()->mesh_vertices.end(),
- v[j]));
- MVertex *v1b;
- if(needParam){
- // parametric coordinates should be set for the vertex ! (this is
- // not 100 % safe yet, so we reserve that operation for high order
- // meshes)
- GPoint gp = gf->closestPoint(SPoint3(v[j]->x(), v[j]->y(), v[j]->z()), guess);
- if(gp.g()){
- v1b = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
- }
- else{
- v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf);
- Msg::Warning("The point was not projected back to the surface (%g %g %g)",
- v[j]->x(), v[j]->y(), v[j]->z());
- }
- }
- else{
- v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf);
- }
-
- gf->mesh_vertices.push_back(v1b);
- numberedV[out.trifacelist[i * 3 + j] - 1] = v1b;
- delete v[j];
- v[j] = v1b;
- }
+ if (out.trifacelist[i * 3 + j] - 1 >= initialSize){
+ printf("aaaaaaaaaaaaaaargh\n");
+ // if(v[j]->onWhat()->dim() == 3){
+ v[j]->onWhat()->mesh_vertices.erase
+ (std::find(v[j]->onWhat()->mesh_vertices.begin(),
+ v[j]->onWhat()->mesh_vertices.end(),
+ v[j]));
+ MVertex *v1b;
+ if(needParam){
+ // parametric coordinates should be set for the vertex ! (this is
+ // not 100 % safe yet, so we reserve that operation for high order
+ // meshes)
+ GPoint gp = gf->closestPoint(SPoint3(v[j]->x(), v[j]->y(), v[j]->z()), guess);
+ if(gp.g()){
+ v1b = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
+ }
+ else{
+ v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf);
+ Msg::Warning("The point was not projected back to the surface (%g %g %g)",
+ v[j]->x(), v[j]->y(), v[j]->z());
+ }
+ }
+ else{
+ v1b = new MVertex(v[j]->x(), v[j]->y(), v[j]->z(), gf);
+ }
+
+ gf->mesh_vertices.push_back(v1b);
+ numberedV[out.trifacelist[i * 3 + j] - 1] = v1b;
+ delete v[j];
+ v[j] = v1b;
+ }
}
// printf("new triangle %d %d %d\n",v[0]->onWhat()->dim(), v[1]->onWhat()->dim(), v[2]->onWhat()->dim());
MTriangle *t = new MTriangle(v[0], v[1], v[2]);
@@ -432,7 +432,7 @@ bool CreateAnEmptyVolumeMesh(GRegion *gr)
buildTetgenStructure(gr, in, numberedV, sqr);
printf("tetgen structure created\n");
sprintf(opts, "-Ype%sT%g",
- (Msg::GetVerbosity() < 3) ? "Q" : (Msg::GetVerbosity() > 6) ? "V" : "",
+ (Msg::GetVerbosity() < 3) ? "Q" : (Msg::GetVerbosity() > 6) ? "V" : "",
CTX::instance()->mesh.toleranceInitialDelaunay);
try{
tetrahedralize(opts, &in, &out);
@@ -952,8 +952,8 @@ void meshGRegion::operator() (GRegion *gr)
if(CTX::instance()->mesh.algo3d == ALGO_3D_FRONTAL){
for(std::list<GFace*>::iterator it = faces.begin(); it != faces.end(); it++){
if((*it)->quadrangles.size()){
- Msg::Error("Cannot use frontal 3D algorithm with quadrangles on boundary");
- return;
+ Msg::Error("Cannot use frontal 3D algorithm with quadrangles on boundary");
+ return;
}
}
}
diff --git a/Mesh/meshGRegionRelocateVertex.cpp b/Mesh/meshGRegionRelocateVertex.cpp
index f913e1033a6a22267ed45b3a53c58a20dbd089d1..32522126bbd8537cd5cb8c0a361e360c42a463ce 100644
--- a/Mesh/meshGRegionRelocateVertex.cpp
+++ b/Mesh/meshGRegionRelocateVertex.cpp
@@ -11,8 +11,8 @@
#include "meshGFaceOptimize.h"
static double objective_function (double xi, MVertex *ver,
- double xTarget, double yTarget, double zTarget,
- const std::vector<MElement*> <){
+ double xTarget, double yTarget, double zTarget,
+ const std::vector<MElement*> <){
double x = ver->x();
double y = ver->y();
double z = ver->z();
@@ -25,6 +25,7 @@ static double objective_function (double xi, MVertex *ver,
minQual = std::min((lt[i]->minSICNShapeMeasure()), minQual);
else
minQual = std::min(fabs(lt[i]->minSICNShapeMeasure()), minQual);
+// minQual = std::min(lt[i]->minAnisotropyMeasure(), minQual);
}
ver->x() = x;
ver->y() = y;
@@ -33,8 +34,8 @@ static double objective_function (double xi, MVertex *ver,
}
static double objective_function (double xi, MVertex *ver, GFace *gf,
- SPoint2 &p1, SPoint2 &p2,
- const std::vector<MElement*> <){
+ SPoint2 &p1, SPoint2 &p2,
+ const std::vector<MElement*> <){
double x = ver->x();
double y = ver->y();
double z = ver->z();
@@ -67,8 +68,8 @@ static int Stopping_Rule(double x0, double x1, double tol)
}
double Maximize_Quality_Golden_Section( MVertex *ver,
- double xTarget, double yTarget, double zTarget,
- const std::vector<MElement*> < , double tol)
+ double xTarget, double yTarget, double zTarget,
+ const std::vector<MElement*> < , double tol)
{
static const double lambda = 0.5 * (sqrt5 - 1.0);
@@ -107,9 +108,9 @@ double Maximize_Quality_Golden_Section( MVertex *ver,
double Maximize_Quality_Golden_Section( MVertex *ver, GFace *gf,
- SPoint2 &p1, SPoint2 &p2,
- const std::vector<MElement*> < ,
- double tol, double &worst)
+ SPoint2 &p1, SPoint2 &p2,
+ const std::vector<MElement*> < ,
+ double tol, double &worst)
{
static const double lambda = 0.5 * (sqrt5 - 1.0);
@@ -154,7 +155,7 @@ double Maximize_Quality_Golden_Section( MVertex *ver, GFace *gf,
}
static void _relocateVertex(MVertex *ver,
- const std::vector<MElement*> <, double tol)
+ const std::vector<MElement*> <, double tol)
{
if(ver->onWhat()->dim() != 3) return;
double x = 0, y=0, z=0;
@@ -179,7 +180,7 @@ static void _relocateVertex(MVertex *ver,
}
static double _relocateVertex(GFace* gf, MVertex *ver,
- const std::vector<MElement*> <, double tol) {
+ const std::vector<MElement*> <, double tol) {
if(ver->onWhat()->dim() != 2) return 2.0;
SPoint2 p1(0,0);
@@ -195,9 +196,9 @@ static double _relocateVertex(GFace* gf, MVertex *ver,
reparamMeshVertexOnFace(v, gf, pp);
counter++;
if (v->onWhat()->dim() == 1) {
- GEdge *ge = dynamic_cast<GEdge*> (v->onWhat());
- // do not take any chance
- if (ge->isSeam(gf))return 2.0;
+ GEdge *ge = dynamic_cast<GEdge*> (v->onWhat());
+ // do not take any chance
+ if (ge->isSeam(gf))return 2.0;
}
p1 += pp;
}