diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 1f9d03dd6265f780c9a43c4a313a153a822a1941..f0d8d20ec5d5c81e093579dc7bd5301e91b52655 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -882,6 +882,74 @@ void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1,
}
+//========================================================================================================
+
+void Curvature::triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs)
+{
+ MVertex* A = triangle->getVertex(0);
+ MVertex* B = triangle->getVertex(1);
+ MVertex* C = triangle->getVertex(2);
+
+ int V0 = 0;
+ int V1 = 0;
+ int V2 = 0;
+
+ std::map<int,int>::iterator vertexIterator;
+ vertexIterator = _VertexToInt.find( A->getNum() );
+ if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
+ else
+ std::cout << "Didn't find vertex with number " << A->getNum() << " in _VertextToInt !" << std::endl;
+
+ vertexIterator = _VertexToInt.find( B->getNum() );
+ if ( vertexIterator != _VertexToInt.end() ) V1 = (*vertexIterator).second;
+ else
+ std::cout << "Didn't find vertex with number " << B->getNum() << " in _VertextToInt !" << std::endl;
+
+ vertexIterator = _VertexToInt.find( C->getNum() );
+ if ( vertexIterator != _VertexToInt.end() ) V2 = (*vertexIterator).second;
+ else
+ std::cout << "Didn't find vertex with number " << C->getNum() << " in _VertextToInt !" << std::endl;
+
+ if (isAbs){
+ dMax[0] = _pdir1[V0];
+ dMax[1] = _pdir1[V1];
+ dMax[2] = _pdir1[V2];
+
+ dMin[0] = _pdir2[V0];
+ dMin[1] = _pdir2[V1];
+ dMin[2] = _pdir2[V2];
+
+ cMax[0] = std::abs(_curv1[V0]);
+ cMax[1] = std::abs(_curv1[V1]);
+ cMax[2] = std::abs(_curv1[V2]);
+
+ cMin[0] = std::abs(_curv2[V0]);
+ cMin[1] = std::abs(_curv2[V1]);
+ cMin[2] = std::abs(_curv2[V2]);
+
+ }
+ else{
+
+ dMax[0] = _pdir1[V0];
+ dMax[1] = _pdir1[V1];
+ dMax[2] = _pdir1[V2];
+
+ dMin[0] = _pdir2[V0];
+ dMin[1] = _pdir2[V1];
+ dMin[2] = _pdir2[V2];
+
+ cMax[0] = _curv1[V0];
+ cMax[1] = _curv1[V1];
+ cMax[2] = _curv1[V2];
+
+ cMin[0] = _curv2[V0];
+ cMin[1] = _curv2[V1];
+ cMin[2] = _curv2[V2];
+ }
+}
+
+
+
//========================================================================================================
void Curvature::edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs)
@@ -922,6 +990,8 @@ void Curvature::writeToPosFile( const std::string & filename)
outfile.open(filename.c_str());
outfile << "View \"Curvature \"{" << std::endl;
+ int idxelem = 0;
+
for (int i = 0; i< _ptFinalEntityList.size(); ++i)
{
GFace* face = _ptFinalEntityList[i]; //face is a pointer to one surface of the group "FinalEntityList"
@@ -944,7 +1014,7 @@ void Curvature::writeToPosFile( const std::string & filename)
//Here is printing the triplet X-Y-Z of each vertex:
//*************************************************
- outfile << "ST("; //VT = vector triangles //ST = scalar triangle
+ outfile << idxelem << " ST("; //VT = vector triangles //ST = scalar triangle
outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
outfile << C->x() << ","<< C->y() << "," << C->z();
@@ -963,7 +1033,10 @@ void Curvature::writeToPosFile( const std::string & filename)
outfile << "};" << std::endl;
+ idxelem++;
+
} //Loop over elements
+ std::cout << "The Total Number of triangle is: " << idxelem << std::endl;
} // Loop over ptFinalEntityList
@@ -1093,7 +1166,114 @@ void Curvature::writeToVtkFile( const std::string & filename)
}
+//========================================================================================================
+
+void Curvature::writeDirectionsToPosFile( const std::string & filename)
+{
+ std::ofstream outfile;
+ outfile.precision(18);
+ outfile.open(filename.c_str());
+ outfile << "View \"Curvature_DirMax \"{" << std::endl;
+
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i)
+ {
+ GFace* face = _ptFinalEntityList[i]; //face is a pointer to one surface of the group "FinalEntityList"
+
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
+ {
+ MElement *e = face->getMeshElement(iElem); //Pointer to one element
+ const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
+
+ //std::cout << "We are now looking at element Nr: " << E << std::endl;
+
+ MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
+ MVertex* B = e->getVertex(1);
+ MVertex* C = e->getVertex(2);
+
+ const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
+ const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
+ const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
+
+ //Here is printing the triplet X-Y-Z of each vertex:
+ //*************************************************
+
+ outfile << "VT("; //VT = vector triangles //ST = scalar triangle
+ outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
+ outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
+ outfile << C->x() << ","<< C->y() << "," << C->z();
+
+ outfile << ")";
+ outfile <<"{";
+
+ //Here is printing the 3 components of the curvature max direction for each vertex:
+ //*********************************************************************************
+
+ outfile << _pdir1[V1].x() << ","<< _pdir1[V1].y() << ","<< _pdir1[V1].z() << ",";
+ outfile << _pdir1[V2].x() << ","<< _pdir1[V2].y() << ","<< _pdir1[V2].z() << ",";
+ outfile << _pdir1[V3].x() << ","<< _pdir1[V3].y() << ","<< _pdir1[V3].z();
+
+
+ outfile << "};" << std::endl;
+
+ } //Loop over elements
+
+} // Loop over ptFinalEntityList
+
+outfile << "};" << std::endl;
+
+
+//----------------------------------------------------------------------------------------------
+
+outfile << "View \"Curvature_DirMin \"{" << std::endl;
+for (int i = 0; i< _ptFinalEntityList.size(); ++i)
+{
+ GFace* face = _ptFinalEntityList[i]; //face is a pointer to one surface of the group "FinalEntityList"
+
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) //Loop over the element all the element of the "myTag"-surface
+ {
+ MElement *e = face->getMeshElement(iElem); //Pointer to one element
+ const int E = _ElementToInt[e->getNum()]; //The NEW tag of the corresponding element
+
+ //std::cout << "We are now looking at element Nr: " << E << std::endl;
+
+ MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
+ MVertex* B = e->getVertex(1);
+ MVertex* C = e->getVertex(2);
+
+ const int V1 = _VertexToInt[A->getNum()]; //Tag of the 1st vertex of the triangle
+ const int V2 = _VertexToInt[B->getNum()]; //Tag of the 2nd vertex of the triangle
+ const int V3 = _VertexToInt[C->getNum()]; //Tag of the 3rd vertex of the triangle
+
+ //Here is printing the triplet X-Y-Z of each vertex:
+ //*************************************************
+
+ outfile << "VT("; //VT = vector triangles //ST = scalar triangle
+ outfile << A->x() << ","<< A->y() << "," << A->z()<< ",";
+ outfile << B->x() << ","<< B->y() << "," << B->z()<< ",";
+ outfile << C->x() << ","<< C->y() << "," << C->z();
+
+ outfile << ")";
+ outfile <<"{";
+
+ //Here is printing the 3 components of the curvature min direction for each vertex:
+ //*********************************************************************************
+
+ outfile << _pdir2[V1].x() << ","<< _pdir2[V1].y() << ","<< _pdir2[V1].z() << ",";
+ outfile << _pdir2[V2].x() << ","<< _pdir2[V2].y() << ","<< _pdir2[V2].z() << ",";
+ outfile << _pdir2[V3].x() << ","<< _pdir2[V3].y() << ","<< _pdir2[V3].z();
+
+
+ outfile << "};" << std::endl;
+
+} //Loop over elements
+
+} // Loop over ptFinalEntityList
+
+outfile << "};" << std::endl;
+
+outfile.close();
+}
//========================================================================================================
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index 898de2be1cb071ce5fe2280c56e584eac873f6f4..60281189580fc2ffb5193fcc4f335eb98a1f6e62 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -178,12 +178,16 @@ public:
void computeCurvature_Rusinkiewicz(int isMax=0);
void triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs=0);
+ void triangleNodalValuesAndDirections(MTriangle* triangle, SVector3* dMax, SVector3* dMin, double* cMax, double* cMin, int isAbs=0);
+
void edgeNodalValues(MLine* edge, double& c0, double& c1, int isAbs=0);
void writeToPosFile( const std::string & filename);
void writeToVtkFile( const std::string & filename);
+ void writeDirectionsToPosFile( const std::string & filename);
+
};
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 03c554194c2c4085a47c51094ef6dc3c4a5d9713..e5f07e70c393b1a5c1620acd44618dbfad05990b 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1417,7 +1417,67 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
return 0.;
}
+double GFaceCompound::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const
+{
+
+ if(!oct) parametrize();
+ if(trivial()) {
+ //Implement this
+// return (*(_compound.begin()))->curvatureMax(param);
+ }
+
+ double U, V;
+ GFaceCompoundTriangle *lt;
+ getTriangle(param.x(), param.y(), <, U,V);
+
+ if(!lt) {
+ return 0.0;
+ }
+
+ if(lt->gf && lt->gf->geomType() != GEntity::DiscreteSurface) {
+ //Implement this...
+// SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
+// return lt->gf->curvatureMax(pv);
+ }
+ else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
+
+ Curvature& curvature = Curvature::getInstance();
+ if( !Curvature::valueAlreadyComputed() ) {
+ std::cout << "Need to compute discrete curvature" << std::endl;
+ std::cout << "Getting instance of curvature" << std::endl;
+
+ curvature.setGModel( model() );
+ int computeMax = 0;
+ curvature.computeCurvature_Rusinkiewicz(computeMax);
+ curvature.writeToPosFile("curvature.pos");
+ curvature.writeToVtkFile("curvature.vtk");
+ curvature.writeDirectionsToPosFile("curvature_directions.pos");
+ std::cout << " ... computing curvature finished" << std::endl;
+ }
+
+ std::cout << "I'm using curvatures in GFaceCompound.cpp" << std::endl;
+ double cMin[3];
+ double cMax[3];
+ SVector3 dMin[3];
+ SVector3 dMax[3];
+
+ curvature.triangleNodalValuesAndDirections(lt->tri, dMax, dMin, cMax, cMin, 0);
+ //curvature.triangleNodalValuesAndDirections(lt->tri, dMax, dMin, cMax, cMin, 1);
+
+ * dirMax = (1-U-V)*dMax[0] + U*dMax[1] + V*dMax[2];
+ * dirMin = (1-U-V)*dMin[0] + U*dMin[1] + V*dMin[2];
+ * curvMax = (1-U-V)*cMax[0] + U*cMax[1] + V*cMax[2];
+ * curvMin = (1-U-V)*cMin[0] + U*cMin[1] + V*cMin[2];
+
+ return * curvMax;
+
+
+ }
+
+ return 0.;
+}
double GFaceCompound::locCurvature(MTriangle *t, double u, double v) const
{
diff --git a/Geo/discreteFace.cpp b/Geo/discreteFace.cpp
index adcf923389b554a864d2ed1e9496d80ed53f32ad..49701ef0a7b58707a8dc1b098260cd1e19b77a9a 100644
--- a/Geo/discreteFace.cpp
+++ b/Geo/discreteFace.cpp
@@ -102,6 +102,18 @@ double discreteFace::curvatureMax(const SPoint2 ¶m) const
}
}
+double discreteFace::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const
+{
+ if (getCompound()){
+ return getCompound()->curvatures(param, dirMax, dirMin, curvMax, curvMin);
+ }
+ else{
+ Msg::Error("Cannot evaluate curvatures and curvature directions on discrete face");
+ return false;
+ }
+}
+
Pair<SVector3, SVector3> discreteFace::firstDer(const SPoint2 ¶m) const
{
if (getCompound()){
diff --git a/Geo/discreteFace.h b/Geo/discreteFace.h
index 446bdbfd10eb3798e5e7b8d7eb2e08a51c488605..6d2980ab578f063486a6805417f0582de3992fb6 100644
--- a/Geo/discreteFace.h
+++ b/Geo/discreteFace.h
@@ -19,6 +19,8 @@ class discreteFace : public GFace {
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
SVector3 normal(const SPoint2 ¶m) const;
double curvatureMax(const SPoint2 ¶m) const;
+ double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
+ double *curvMax, double *curvMin) const;
GEntity::GeomType geomType() const { return DiscreteSurface; }
virtual Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 ¶m,