diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp index 3bfa04943ec4bb4a2d57e8f967f5997e86fe235d..cc2a5d2d196d55a0e9e354845c853351a0f2c4ac 100644 --- a/Mesh/directions3D.cpp +++ b/Mesh/directions3D.cpp @@ -4,7 +4,7 @@ // bugs and problems to the public mailing list <gmsh@geuz.org>. // // Contributor(s): -// Tristan Carrier +// Tristan Carrier François Henrotte #include <fstream> #include "GModel.h" @@ -95,23 +95,23 @@ void Frame_field::init_face(GFace* gf){ for(j=0;j<element->getNumVertices();j++){ vertex = element->getVertex(j); - ok = translate(gf,octree,vertex,SPoint2(average_x,average_y),v1,v2); + ok = improved_translate(gf,vertex,v1,v2); if(ok){ - 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()); - temp.insert(std::pair<MVertex*,Matrix>(vertex,m)); + 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()); + temp.insert(std::pair<MVertex*,Matrix>(vertex,m)); } } } @@ -173,6 +173,31 @@ bool Frame_field::translate(GFace* gf,MElementOctree* octree,MVertex* vertex,SPo return ok; } +bool Frame_field::improved_translate(GFace* gf,MVertex* vertex,SVector3& v1,SVector3& v2){ + double x,y; + double angle; + SPoint2 point; + SVector3 s1,s2; + SVector3 normal; + 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); + + v1 = s1*cos(angle) + s2*sin(angle); + v2 = crossprod(v1,normal); + + return 1; +} + Matrix Frame_field::search(double x,double y,double z){ // Determines the frame/cross at location (x,y,z) int index=0; @@ -807,8 +832,9 @@ void Size_field::solve(GRegion* gr){ //printf("number of tetrahedra = %d\n",count2); //printf("volume = %f\n",volume); - if (assembler.sizeOfR()) + if(assembler.sizeOfR()){ system->systemSolve(); + } for(it=interior.begin();it!=interior.end();it++){ assembler.getDofValue(*it,0,1,val);