diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index 91800deeadef74ab628ae67f6b3e3f2325c953f4..81e03d992f15fdd55a18a16174cec322f071ddc4 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -212,10 +212,8 @@ void statisticsWindow::compute(bool elementQuality)
GModel::current()->getEntities(entities);
std::map<MVertex*, int > vert2Deg;
for(unsigned int i = 0; i < entities.size(); i++){
- printf("entity dim =%d tag=%d \n", entities[i]->dim() , entities[i]->tag());
if(entities[i]->dim() < 2 ) continue;
if(entities[i]->tag() != 10) continue;
- printf("continuing \n");
for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
MElement *e = entities[i]->getMeshElement(j);
for(unsigned int k = 0; k < e->getNumEdges(); k++){
@@ -251,15 +249,13 @@ void statisticsWindow::compute(bool elementQuality)
FieldManager *fields = GModel::current()->getFields();
Field *f = fields->get(fields->background_field);
int nbEdges = edges.size();
- printf("nb edges =%d \n", nbEdges);
system("rm qualEdges.txt");
FILE *fp = fopen("qualEdges.txt", "w");
std::vector<int> qualE;
int nbS = 50;
qualE.resize(nbS);
if(fields->background_field > 0){
- printf("found field \n");
- std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+ std::set<MEdge, Less_Edge>::iterator it = edges.begin();
double sum = 0;
for (; it !=edges.end();++it){
MVertex *v0 = it->getVertex(0);
@@ -273,7 +269,6 @@ void statisticsWindow::compute(bool elementQuality)
sum += e - 1.0;
}
double tau = exp ((1./edges.size()) * sum);
- printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
double ibegin = 2./(2*nbS);
double inext = 2./nbS;
for (int i= 0; i< qualE.size(); i++){
diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index a7538e739965ad7b8fa8fa1ff40a7eb90ced0c8e..05195fbba0b7bd22bce12cfdf11b7ceb9aba92c4 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -10,6 +10,7 @@
#include "GFaceCompound.h"
#include "MLine.h"
#include "GRbf.h"
+#include "Os.h"
#include "SBoundingBox3d.h"
#include<iostream>
@@ -81,13 +82,6 @@ Curvature& Curvature::getInstance()
_instance = & instance;
}
-//========================================================================================================
-
- void Curvature::setGModel(GModel* model)
- {
- _model = model;
- }
-
//========================================================================================================
void Curvature::retrieveCompounds() {
#if defined(HAVE_SOLVER)
@@ -415,10 +409,8 @@ void Curvature::computeRusinkiewiczNormals()
{
// Pointer to one element
MElement *e = face->getMeshElement(iElem);
- // The NEW tag of the corresponding element
const int E = _ElementToInt[e->getNum()];
- // std::cout << "We are now looking at element Nr: " << E << std::endl;
-
+
// Pointers to vertices of triangle
MVertex* A = e->getVertex(0);
MVertex* B = e->getVertex(1);
@@ -432,13 +424,10 @@ void Curvature::computeRusinkiewiczNormals()
vector_AC = SVector3(C->x() - A->x(), C->y() - A->y(), C->z() - A->z() );
vector_BC = SVector3(C->x() - B->x(), C->y() - B->y(), C->z() - B->z() );
-
const SVector3 cross = crossprod(vector_AB, vector_AC);
// Area of the triangles:
_TriangleArea[E] = 0.5*cross.norm();
- // std::cout << "The area of the triangle nr: " << e->getNum() << " is: "<< TriangleArea[E] << std::endl;
-
const double l_AB = vector_AB.normSq();
const double l_AC = vector_AC.normSq();
@@ -460,13 +449,8 @@ void Curvature::computeRusinkiewiczNormals()
}
//========================================================================================================
-
// Compute per-vertex point areas
-void Curvature::computePointareas()
-{
-
-// std::cout << "Computing point areas... " << std::endl;
-// std::cout << "The mesh has " << _VertexToInt.size() << " nodes" << std::endl;
+void Curvature::computePointareas(){
SVector3 e[3];
SVector3 l2;
@@ -547,17 +531,13 @@ void Curvature::computePointareas()
} //End of loop over _ptFinalEntityList
- //std::cout << "Done computing pointareas" << std::endl;
-
} //End of the method "computePointareas"
//========================================================================================================
-
//Rotate a coordinate system to be perpendicular to the given normal
+void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v, const SVector3 &new_norm, SVector3 &new_u, SVector3 &new_v){
-void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v, const SVector3 &new_norm, SVector3 &new_u, SVector3 &new_v)
-{
new_u = old_u;
new_v = old_v;
SVector3 old_norm = crossprod(old_u, old_v);
@@ -585,8 +565,7 @@ void Curvature::rot_coord_sys(const SVector3 &old_u, const SVector3 &old_v, cons
void Curvature::proj_curv( const SVector3 &old_u, const SVector3 &old_v,
double old_ku, double old_kuv, double old_kv,
const SVector3 &new_u, const SVector3 &new_v,
- double &new_ku, double &new_kuv, double &new_kv)
-{
+ double &new_ku, double &new_kuv, double &new_kv){
SVector3 r_new_u;
SVector3 r_new_v;
rot_coord_sys(new_u, new_v, crossprod(old_u,old_v), r_new_u, r_new_v);
@@ -610,8 +589,7 @@ void Curvature::proj_curv( const SVector3 &old_u, const SVector3 &old_v,
void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
double ku, double kuv, double kv,
const SVector3 &new_norm,
- SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2)
-{
+ SVector3 &pdir1, SVector3 &pdir2, double &k1, double &k2){
SVector3 r_old_u;
SVector3 r_old_v;
@@ -646,8 +624,30 @@ void Curvature::diagonalize_curv(const SVector3 &old_u, const SVector3 &old_v,
}
//========================================================================================================
+void Curvature::computeCurvature(GModel* model, typeOfCurvature typ)
+{
+
+ _model = model;
+
+ double t0 = Cpu();
+ Msg::StatusBar(2, true, "(C) Computing Curvature");
+ if (typ == RUSIN)
+ computeCurvature_Rusinkiewicz(0);
+ else if (typ == RBF)
+ computeCurvature_RBF();
+ else if (typ == SIMPLE)
+ computeCurvature_Simple();
+
+ double t1 = Cpu();
+ Msg::StatusBar(2, true, "(C) Done Computing Curvature (%g s)", t1-t0);
-void Curvature::computeCurvature_Rusinkiewicz(int isMax){
+ writeToPosFile("curvature.pos");
+ writeToVtkFile("curvature.vtk");
+
+}
+
+void Curvature::computeCurvature_Rusinkiewicz(int isMax)
+{
retrieveCompounds();
initializeMap();
computeRusinkiewiczNormals();
@@ -679,8 +679,6 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax){
_curv2.resize(_VertexToInt.size());
_curv12.resize(_VertexToInt.size());
- std::cout << "Computing Curvature.." << 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"
@@ -813,18 +811,14 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax){
//Compute principal directions and curvatures at each vertex
- for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
- {
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex) {
diagonalize_curv(_pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex], _curv12[ivertex], _curv2[ivertex],
_VertexNormal[ivertex], _pdir1[ivertex], _pdir2[ivertex], _curv1[ivertex], _curv2[ivertex]);
}
- std::cout << "Done" << std::endl;
-
_VertexCurve.resize( _VertexToInt.size() );
- for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex)
- {
+ for (int ivertex = 0; ivertex < _VertexToInt.size(); ++ivertex){
if (isMax){
_VertexCurve[ivertex] = std::max(fabs(_curv1[ivertex]), fabs(_curv2[ivertex]));
@@ -842,11 +836,10 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax){
} //End of the "computeCurvature_Rusinkiewicz" method
-void Curvature::computeCurvature_RBF(){
+void Curvature::computeCurvature_RBF()
+{
retrieveCompounds();
initializeMap();
-
- std::cout << "Computing Curvature RBF ..." << std::endl;
//fill set of MVertex
std::set<MVertex*> allNodes;
@@ -869,17 +862,17 @@ void Curvature::computeCurvature_RBF(){
bb +=pt;
}
double sizeBox = norm(SVector3(bb.max(), bb.min()));
- printf("allNodes := %d sizeBox = %g \n", allNodes.size(), sizeBox);
//compure curvature RBF
std::map<MVertex*, SVector3> _normals;
std::vector<MVertex*> _ordered;
std::map<MVertex*, double> curvRBF;
- GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered); //, true);
+ //GLOBAL
+ GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered);
_rbf->computeCurvature(_rbf->getXYZ(),curvRBF);
- // _rbf->computeLocalCurvature(_rbf->getXYZ(),curvRBF);
-
- std::cout << "Done" << std::endl;
+ //LOCAL FD
+ //GRbf *_rbf = new GRbf(sizeBox, 0, 1, _normals, allNodes, _ordered, true);
+ //_rbf->computeLocalCurvature(_rbf->getXYZ(),curvRBF);
//fill vertex curve
_VertexCurve.resize( _VertexToInt.size() );
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index 349626ac34dd1d201e8c8f342b459f1b3efdd800..206f370cef15db9a9d7fd05b38c1e3c220c955f8 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -13,11 +13,8 @@
#include<vector>
class Curvature {
-private:
-
- //-----------------------------------------
- // TYPEDEFS
+private:
typedef std::vector<GFace*> GFaceList;
//typedef std::map<int, GEntityList >::iterator GEntityIter;
typedef std::map<int, GFaceList >::iterator GFaceIter;
@@ -102,7 +99,6 @@ private:
// Perform LDL^T decomposition of a symmetric positive definite matrix.
// Like Cholesky, but no square roots. Overwrites lower triangle of matrix.
-
static inline bool ldltdc(STensor3& A, double rdiag[3])
{
double v[2];
@@ -131,7 +127,6 @@ private:
return true;
}
-
// Solve Ax=B after ldltdc
static inline void ldltsl(STensor3& A, double rdiag[3], double B[3], double x[3])
{
@@ -152,20 +147,18 @@ private:
}
}
-
public:
- //-----------------------------------------
- // PUBLIC METHODS
-
+ typedef enum {RUSIN=1,RBF=2, SIMPLE=3} typeOfCurvature;
static Curvature& getInstance();
static bool valueAlreadyComputed();
- void setGModel(GModel* model);
void retrieveCompounds();
double getAtVertex(const MVertex *v) const;
//void retrievePhysicalSurfaces(const std::string & face_tag);
+ void computeCurvature(GModel* model, typeOfCurvature typ);
+
/// The following function implements algorithm from:
/// Implementation of an Algorithm for Approximating the Curvature Tensor
/// on a Triangular Surface Mesh in the Vish Environment
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index c2f7e79a8fca51301a7da7c80dab6f6332386968..f5c0b72b1439371c619632155e552406e1b53dd1 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -658,9 +658,9 @@ bool GFaceCompound::parametrize() const
_rbf = new GRbf(sizeBox, variableEps, radFunInd, _normals, allNodes, _ordered);
//_rbf->RbfLapSurface_global_CPM_low(_rbf->getXYZ(), _rbf->getN(), Oper);
- //_rbf->RbfLapSurface_local_CPM(true, _rbf->getXYZ(), _rbf->getN(), Oper, true);
- _rbf->RbfLapSurface_global_CPM_high(_rbf->getXYZ(), _rbf->getN(), Oper);
- //_rbf->RbfLapSurface_local_CPM(false, _rbf->getXYZ(), _rbf->getN(), Oper, true);
+ //_rbf->RbfLapSurface_local_CPM(true, _rbf->getXYZ(), _rbf->getN(), Oper);
+ _rbf->RbfLapSurface_global_CPM_high_2(_rbf->getXYZ(), _rbf->getN(), Oper);
+ //_rbf->RbfLapSurface_local_CPM(false, _rbf->getXYZ(), _rbf->getN(), Oper);
//_rbf->RbfLapSurface_global_projection(_rbf->getXYZ(), _rbf->getN(), Oper);
//_rbf->RbfLapSurface_local_projection(_rbf->getXYZ(), _rbf->getN(), Oper, true);
@@ -1379,30 +1379,16 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
return lt->gf->curvatureMax(pv);
}
else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
-
Curvature& curvature = Curvature::getInstance();
-
if( !Curvature::valueAlreadyComputed() ) {
- Msg::Info("Need to compute discrete curvature for isotropic remesh");
- Msg::Info("Getting instance of curvature!");
-
- curvature.setGModel( model() );
- int computeMax = 0;
- curvature.computeCurvature_Rusinkiewicz(computeMax);
- //curvature.computeCurvature_RBF();
- curvature.writeToPosFile("curvature.pos");
- curvature.writeToVtkFile("curvature.vtk");
- Msg::Info(" ... computing curvature finished");
+ Msg::Info("Need to compute discrete curvature for isotropic remesh (in GFace)");
+ Curvature::typeOfCurvature type = Curvature::RUSIN;
+ curvature.computeCurvature(model(), type);
}
-
- double c0;
- double c1;
- double c2;
+ double c0,c1,c2;
curvature.triangleNodalValues(lt->tri,c0, c1, c2, 1);
-
double cv = (1-U-V)*c0 + U*c1 + V*c2;
return cv;
-
}
return 0.;
@@ -1413,8 +1399,7 @@ double GFaceCompound::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector
if(!oct) parametrize();
if(trivial()) {
- //Implement this
-// return (*(_compound.begin()))->curvatureMax(param);
+ return (*(_compound.begin()))->curvatures(param, dirMax,dirMin, curvMax, curvMin);
}
double U, V;
@@ -1426,34 +1411,22 @@ double GFaceCompound::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector
}
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);
+ SPoint2 pv = lt->gfp1*(1.-U-V) + lt->gfp2*U + lt->gfp3*V;
+ return lt->gf->curvatures(pv, dirMax,dirMin, curvMax, curvMin);
}
- else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
+ else if (lt->gf->geomType() == GEntity::DiscreteSurface) {
Curvature& curvature = Curvature::getInstance();
-
if( !Curvature::valueAlreadyComputed() ) {
- Msg::Info("Need to compute discrete curvature for anisotropic remesh");
- Msg::Info("Getting instance of curvature");
-
- curvature.setGModel( model() );
- int computeMax = 0;
- curvature.computeCurvature_Rusinkiewicz(computeMax);
- //curvature.computeCurvature_RBF();
- curvature.writeToPosFile("curvature.pos");
- //curvature.writeToVtkFile("curvature.vtk");
- //curvature.writeDirectionsToPosFile("curvature_directions.pos");
- Msg::Info(" ... computing curvature finished");
+ Msg::Info("Need to compute discrete curvature for anisotropic remesh (in GFace)");
+ Curvature::typeOfCurvature type = Curvature::RUSIN;//RBF
+ curvature.computeCurvature(model(), type);
}
- 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);
@@ -1463,8 +1436,6 @@ double GFaceCompound::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector
* curvMin = (1-U-V)*cMin[0] + U*cMin[1] + V*cMin[2];
return * curvMax;
-
-
}
return 0.;
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index e921d4a7d79fef4582a2aa61d2d22d182a30375c..cd8af08eff2cce8b1a0728f08877f4c39eb37ef2 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -78,7 +78,7 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
#endif
allCenters.resize(allNodes.size(),3);
- double tol = 6.e-1*sBox;
+ double tol = 4.e-2*sBox;
if (isLocal) tol = 1.e-15;
//remove duplicate vertices
@@ -203,13 +203,13 @@ void GRbf::buildOctree(double radius){
double P[3] = {centers(i,0),centers(i,1), centers(i,2)};
Octree_SearchAll(P, oct, &l);
nodesInSphere[i].push_back(i);
- if (l.size()== 1) printf("*** WARNING: Found only %d sphere ! \n", (int)l.size());
+ if (l.size() == 1) printf("*** WARNING: Found only %d sphere ! \n", (int)l.size());
for (std::list<void*>::iterator it = l.begin(); it != l.end(); it++) {
Sphere *sph = (Sphere*) *it;
std::vector<int> &pts = nodesInSphere[i];
if (sph->index != i) nodesInSphere[i].push_back(sph->index);
}
- printf("size node i =%d = %d \n", i , nodesInSphere[i].size());
+ //printf("size node i =%d = %d \n", i , nodesInSphere[i].size());
}
Octree_Delete(oct);
@@ -228,18 +228,12 @@ void GRbf::curvatureRBF(const fullMatrix<double> &cntrs,
evalRbfDer(1,extX,cntrs,surf,sx);
evalRbfDer(2,extX,cntrs,surf,sy);
evalRbfDer(3,extX,cntrs,surf,sz);
- //evalRbfDer(11,extX,cntrs,surf,sxx);
- //evalRbfDer(12,extX,cntrs,surf,sxy);
- //evalRbfDer(13,extX,cntrs,surf,sxz);
- //evalRbfDer(22,extX,cntrs,surf,syy);
- //evalRbfDer(23,extX,cntrs,surf,syz);
- //evalRbfDer(33,extX,cntrs,surf,szz);
evalRbfDer(222,extX,cntrs,surf,sLap);
for (int i = 0; i < cntrs.size1(); i++) {
double norm_grad_s = sqrt(sx(i,0)*sx(i,0)+sy(i,0)*sy(i,0)+sz(i,0)*sz(i,0));
- double curv = -sLap(i,0)/norm_grad_s; //+(sx(i,0)*sx(i,0)*sxx(i,0)+sy(i,0)*sy(i,0)*syy(i,0)+sz(i,0)*sz(i,0)*szz(i,0)+2*sx(i,0)*sy(i,0)*sxy(i,0)+2*sy(i,0)*sz(i,0)*syz(i,0)+2*sx(i,0)*sz(i,0)*sxz(i,0))/ (norm_grad_s*norm_grad_s*norm_grad_s);
- curvature(i,0) = fabs(curv)/sBox;
+ double curv = -sLap(i,0)/norm_grad_s;
+ curvature(i,0) = 0.5*fabs(curv)/sBox;
}
}
@@ -284,7 +278,6 @@ void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
fullMatrix<double> curv(pts.size(), 1);
curvatureRBF(nodes_in_sph,curv);
- printf("curv(0,0) = %g \n", curv(0,0));
curvature(i,0) = curv(0,0);
}
@@ -339,7 +332,7 @@ fullMatrix<double> GRbf::generateRbfMat(int p,
int n = nodes1.size1();
fullMatrix<double> rbfMat(m,n);
- double eps =0.5/delta; //0.0677*(nbNodes^0.28)/min_dist; //0.5
+ double eps = 0.5/delta; //0.0677*(nbNodes^0.28)/min_dist; //0.5
if (_inUV) eps = 0.4/deltaUV;
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
@@ -411,12 +404,12 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
//Computes the normal vectors to the surface at each node
for (int i=0;i<numNodes ; ++i){
- ONES(i,0)=1.0;
+ ONES(i,0) = 0.0;
cntrsPlus(i,0) = cntrs(i,0);
cntrsPlus(i,1) = cntrs(i,1);
cntrsPlus(i,2) = cntrs(i,2);
}
- ONES(numNodes,0) = 5.0;
+ ONES(numNodes,0) = 1.0;
cntrsPlus(numNodes,0) = cntrs(0,0)+10.*sBox;
cntrsPlus(numNodes,1) = cntrs(0,1)+10.*sBox;
cntrsPlus(numNodes,2) = cntrs(0,2)+10.*sBox;
@@ -425,6 +418,11 @@ void GRbf::setup_level_set(const fullMatrix<double> &cntrs,
evalRbfDer(2,cntrsPlus,cntrs,ONES,sy);
evalRbfDer(3,cntrsPlus,cntrs,ONES,sz);
for (int i=0;i<numNodes ; ++i){
+ //GMSH NORMALS
+ //sx(i,0) = normals(i,0);
+ //sy(i,0) = normals(i,1);
+ //sz(i,0) = normals(i,2);
+ //END GMSH NORMALS
normFactor = sqrt(sx(i,0)*sx(i,0)+sy(i,0)*sy(i,0)+sz(i,0)*sz(i,0));
sx(i,0) = sx(i,0)/normFactor;
sy(i,0) = sy(i,0)/normFactor;
@@ -562,7 +560,7 @@ void GRbf::RbfLapSurface_local_CPM(bool isLow,
LocalOper.setAll(0.0);
if (isLow) RbfLapSurface_global_CPM_low(nodes_in_sph,local_normals,LocalOper);
- else RbfLapSurface_global_CPM_high(nodes_in_sph,local_normals,LocalOper);
+ else RbfLapSurface_global_CPM_high_2(nodes_in_sph,local_normals,LocalOper);
for (int j = 0; j < nbp ; j++){
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index c46aed378dcac3a5ecb236714f98f6839ff0db8d..59948c017966a9ddb3cec68dfdfb356cfa9dd852 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -499,22 +499,11 @@ double discreteEdge::curvature(double par) const{
if(!getLocalParameter(par,iEdge,tLoc)) return MAX_LC;
double c0, c1;
-
Curvature& curvature = Curvature::getInstance();
-
if( !Curvature::valueAlreadyComputed() ) {
- std::cout << "Need to compute discrete curvature (in Edge)" << std::endl;
- std::cout << "Getting instance of curvature" << std::endl;
-
- curvature.setGModel( model() );
- int computeMax = 0;
- curvature.computeCurvature_Rusinkiewicz(computeMax);
- //curvature.computeCurvature_RBF();
- curvature.writeToPosFile("curvature.pos");
- //curvature.writeDirectionsToPosFile("curvature_directions.pos");
- //curvature.writeToVtkFile("curvature.vtk");
-
- std::cout << " ... computing curvature finished" << std::endl;
+ std::cout << "Need to compute discrete curvature (in discreteEdge)" << std::endl;
+ Curvature::typeOfCurvature type = Curvature::RUSIN; //RBF
+ curvature.computeCurvature(model(), type);
}
curvature.edgeNodalValues(lines[iEdge],c0, c1, 1);
@@ -525,7 +514,6 @@ double discreteEdge::curvature(double par) const{
}
-
Range<double> discreteEdge::parBounds(int i) const
{
return Range<double>(0, lines.size());
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 1edcd97bd388b128ab81abdbad6774b4abaef8fe..2f18c7c87e344a66903e56e5e3fe8f8a25e436fb 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -482,8 +482,8 @@ class functionLevelsetSmooth : public function {
//double Heps = 0.5 * (1 + phi / _E + 1. / M_PI * sin(M_PI * phi / _E));
//double Heps = 0.5 + 1./32.*(45.*phi/_E - 50.*pow(phi/_E,3.) + 21.*pow(phi/_E,5.) );
- //double Heps = 0.5*(1+tanh(M_PI*phi/_E));
- double Heps = 0.75 * (phi/_E - 0.33*pow(phi/_E,3.0)) + 0.5;
+ double Heps = 0.5*(1+tanh(M_PI*phi/_E));
+ //double Heps = 0.75 * (phi/_E - 0.33*pow(phi/_E,3.0)) + 0.5;
//if (fabs(phi) < _E) val(i, j) = 1./(Heps * ivalPlus + (1 - Heps) * ivalMin);
//else if (phi > _E) val(i, j) = 1./ivalPlus;
diff --git a/benchmarks/extrude/aorta_extrude.geo b/benchmarks/extrude/aorta_extrude.geo
index e4197d7384d1a302f7890c22f716ffa59eaf5fd6..58b7881b7da531d3aaba23d510224e7fbe253c75 100644
--- a/benchmarks/extrude/aorta_extrude.geo
+++ b/benchmarks/extrude/aorta_extrude.geo
@@ -4,7 +4,7 @@ CreateTopology;
//Merge "aortaRADIUS2.bgm";
// create a boundary layer, whose tickness is given in View[0]
-out1[] = Extrude{Surface{1}; Layers{4, 0.5}; Using Index[0]; Using View[0]; };
+out1[] = Extrude{Surface{-1}; Layers{4, 0.5}; Using Index[0]; Using View[0]; };
// we could create a second boundary layer inside...
//out2[] = Extrude{Surface{1}; Layers{4, -0.5}; Using Index[1]; Using View[0]; };
diff --git a/benchmarks/extrude/u_shape_boundary_layer.geo b/benchmarks/extrude/u_shape_boundary_layer.geo
index c256c0f452b7a172813d5ce3a38a96f756bf460a..cb251e96e30d91694c1b264eb213bfaf85edde03 100644
--- a/benchmarks/extrude/u_shape_boundary_layer.geo
+++ b/benchmarks/extrude/u_shape_boundary_layer.geo
@@ -1,3 +1,5 @@
+Mesh.Algorithm = 5; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg)
+
Point(1) = {0, 0.5, 0, 0.1};
Point(2) = {0.1, 0.7, 0, 0.1};
Point(3) = {0.3, 0.8, 0, 0.1};
@@ -11,10 +13,10 @@ BSpline(1) = {5, 4, 3, 2, 1};
BSpline(2) = {1, 9, 8, 7, 5};
Line(3) = {5, 6};
-Extrude { Line{1,-3}; Layers{5,0.04}; Using Index[0]; }
-Extrude { Line{2,3}; Layers{5,0.04}; Using Index[1]; }
+Extrude { Line{1,-3}; Layers{2,0.04}; Using Index[0]; }
+Extrude { Line{2,3}; Layers{2,0.04}; Using Index[1]; }
-// fix leading edge by hand
+//fix leading edge by hand
Coherence Point {25, 16};
Point(31) = {-0.5, 1.5, 0, 0.2};
diff --git a/tutorial/t11.geo b/tutorial/t11.geo
index 68e36af23743a46fbff109a59512a8caf0797fbd..a113e50958d8d0e6b7fcbfdaea516fd4bbdf4b13 100644
--- a/tutorial/t11.geo
+++ b/tutorial/t11.geo
@@ -45,4 +45,4 @@ Recombine Surface{100};
// triangulation algorithm that enables to create right triangles
// almost everywhere. Uncomment the following line to try DelQuad:
-// Mesh.Algorithm = 8; // DelQuad (experimental)
+Mesh.Algorithm = 8; // DelQuad (experimental)
diff --git a/tutorial/t12.geo b/tutorial/t12.geo
index 1021698ee05aa2efc55676c8fff69925ecaf6def..3cf793aead82b3a90fd284be0416783f5c555186 100644
--- a/tutorial/t12.geo
+++ b/tutorial/t12.geo
@@ -38,7 +38,7 @@ Compound Surface(200) = {12, 14, 16};
// Hide the original surfaces so we only see the compound
// (cross-patch) mesh
-Hide {Surface{12, 14, 16}; }
+//Hide {Surface{12, 14, 16}; }
// More details about the reparametrization technique can be found in
// the following papers: