From 1aaad4fab5ed0c1c87bd8e98c99b71a9bde26243 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Tue, 13 Sep 2011 13:37:32 +0000
Subject: [PATCH] multDXiDXJRight unrolled for dPsiPsi fixed ! bubble
benchmark for level set !
---
Geo/Curvature.cpp | 82 ++++++++++---
Geo/Curvature.h | 2 +
Geo/GFaceCompound.cpp | 19 +--
Geo/GRbf.cpp | 114 ++++++++++++------
Geo/GRbf.h | 12 +-
Geo/discreteEdge.cpp | 7 +-
benchmarks/curvature/test.geo | 3 +-
benchmarks/stl/reparam.geo | 6 +-
.../{reparam_input.msh => reparamINPUT.msh} | 0
9 files changed, 175 insertions(+), 70 deletions(-)
rename benchmarks/stl/{reparam_input.msh => reparamINPUT.msh} (100%)
diff --git a/Geo/Curvature.cpp b/Geo/Curvature.cpp
index 83689d1e3f..a7538e7399 100644
--- a/Geo/Curvature.cpp
+++ b/Geo/Curvature.cpp
@@ -9,6 +9,8 @@
#include "GEntity.h"
#include "GFaceCompound.h"
#include "MLine.h"
+#include "GRbf.h"
+#include "SBoundingBox3d.h"
#include<iostream>
#include<fstream>
@@ -121,20 +123,12 @@ void Curvature::initializeMap()
for (int i = 0; i< _ptFinalEntityList.size(); ++i)
{
- // face is a pointer to one surface of the group "FinalEntityList"
GFace* face = _ptFinalEntityList[i];
-// std::cout << "Face " << i << " has " << face->getNumMeshElements() << " elements" << std::endl;
-
- // Loop over the element all the element of the "myTag"-surface
for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++)
{
-
- // Pointer to one element
MElement *e = face->getMeshElement(iElem);
- // Tag of the corresponding element
const int E = e->getNum();
- // std::cout << "We are now looking at element Nr: " << E << std::endl;
_ElementToInt[E] = 1;
const int A = e->getVertex(0)->getNum(); //Pointer to 1st vertex of the triangle
@@ -847,6 +841,62 @@ void Curvature::computeCurvature_Rusinkiewicz(int isMax){
} //End of the "computeCurvature_Rusinkiewicz" method
+
+void Curvature::computeCurvature_RBF(){
+ retrieveCompounds();
+ initializeMap();
+
+ std::cout << "Computing Curvature RBF ..." << std::endl;
+
+ //fill set of MVertex
+ std::set<MVertex*> allNodes;
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i) {
+ GFaceCompound* face = (GFaceCompound*)_ptFinalEntityList[i];
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++) {
+ MElement *e = face->getMeshElement(iElem);
+ for(int j = 0; j < e->getNumVertices(); j++){
+ allNodes.insert(e->getVertex(j));
+ }
+ }
+ }
+
+ //bounding box
+ SBoundingBox3d bb;
+ std::vector<SPoint3> vertices;
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ SPoint3 pt((*itv)->x(),(*itv)->y(), (*itv)->z());
+ vertices.push_back(pt);
+ 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);
+ _rbf->computeCurvature(_rbf->getXYZ(),curvRBF);
+ // _rbf->computeLocalCurvature(_rbf->getXYZ(),curvRBF);
+
+ std::cout << "Done" << std::endl;
+
+ //fill vertex curve
+ _VertexCurve.resize( _VertexToInt.size() );
+ for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ std::map<int,int>::iterator vertexIterator;
+ int V0;
+ vertexIterator = _VertexToInt.find( v->getNum() );
+ if ( vertexIterator != _VertexToInt.end() ) V0 = (*vertexIterator).second;
+ _VertexCurve[V0] = curvRBF[v];
+ }
+
+ _alreadyComputedCurvature = true;
+
+} //End of the "computeCurvature_RBF" method
+
+
//========================================================================================================
void Curvature::triangleNodalValues(MTriangle* triangle, double& c0, double& c1, double& c2, int isAbs)
@@ -1007,15 +1057,12 @@ void Curvature::writeToPosFile( const std::string & filename)
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"
-
- 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
+ for (int i = 0; i< _ptFinalEntityList.size(); ++i) {
+ GFace* face = _ptFinalEntityList[i];
+ for (int iElem = 0; iElem < face->getNumMeshElements(); iElem++){
+ MElement *e = face->getMeshElement(iElem);
+ const int E = _ElementToInt[e->getNum()];
//std::cout << "We are now looking at element Nr: " << E << std::endl;
MVertex* A = e->getVertex(0); //Pointers to vertices of triangle
@@ -1051,8 +1098,7 @@ void Curvature::writeToPosFile( const std::string & filename)
idxelem++;
} //Loop over elements
- std::cout << "The Total Number of triangle is: " << idxelem << std::endl;
-
+
} // Loop over ptFinalEntityList
outfile << "};" << std::endl;
diff --git a/Geo/Curvature.h b/Geo/Curvature.h
index 07a024a412..349626ac34 100644
--- a/Geo/Curvature.h
+++ b/Geo/Curvature.h
@@ -178,6 +178,8 @@ public:
/// Code taken from Rusinkiewicz' 'trimesh2' library
void computeCurvature_Rusinkiewicz(int isMax=0);
+ void computeCurvature_RBF();
+
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);
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index fc00182526..935ff58135 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -612,7 +612,6 @@ bool GFaceCompound::parametrize() const
if(allNodes.empty()) buildAllNodes();
bool success = orderVertices(_U0, _ordered, _coords);
- printBound(_ordered, this->tag());
if(!success) {Msg::Error("Could not order vertices on boundary");exit(1);}
// Laplace parametrization
@@ -659,17 +658,18 @@ 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);
+ //_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);
+ //_rbf->RbfLapSurface_local_CPM(false, _rbf->getXYZ(), _rbf->getN(), Oper, true);
//_rbf->RbfLapSurface_global_projection(_rbf->getXYZ(), _rbf->getN(), Oper);
- //_rbf->RbfLapSurface_local_projection(_rbf->getXYZ(), _rbf->getN(), Oper);
+ //_rbf->RbfLapSurface_local_projection(_rbf->getXYZ(), _rbf->getN(), Oper, true);
_rbf->solveHarmonicMap(Oper, _ordered, _coords, coordinates);
//_rbf->computeCurvature(coordinates);
//printStuff();
//exit(1);
+
double t2 = Cpu();
printf("param RBF in %g sec \n", t2-t1);
}
@@ -1144,6 +1144,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
Msg::Error("Switch to harmonic map or see doc on the wiki for installing petsc and slepc:");
Msg::Error("https://geuz.org/trac/gmsh/wiki/STLRemeshing (username:gmsh,passwd:gmsh)");
return false;
+ parametrize_conformal();
#else
linearSystem <double> *lsysA = new linearSystemPETSc<double>;
linearSystem <double> *lsysB = new linearSystemPETSc<double>;
@@ -1383,11 +1384,12 @@ double GFaceCompound::curvatureMax(const SPoint2 ¶m) const
if( !Curvature::valueAlreadyComputed() ) {
Msg::Info("Need to compute discrete curvature for isotropic remesh");
- Msg::Info("Getting instance of curvature");
+ 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");
@@ -1439,13 +1441,14 @@ double GFaceCompound::curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector
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");
+ //curvature.writeToVtkFile("curvature.vtk");
+ //curvature.writeDirectionsToPosFile("curvature_directions.pos");
Msg::Info(" ... computing curvature finished");
}
- //std::cout << "I'm using curvatures in GFaceCompound.cpp" << std::endl;
+ std::cout << "I'm using curvatures in GFaceCompound.cpp" << std::endl;
double cMin[3];
double cMax[3];
SVector3 dMin[3];
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index 52f9610d75..ffaa761cd9 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -69,12 +69,13 @@ static void exportParametrizedMesh(fullMatrix<double> &UV, int nbNodes){
}
GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVector3> _normals,
- std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes)
- : sBox(sizeBox), variableShapeParam(variableEps), radialFunctionIndex (rbfFun), XYZkdtree(0), _inUV(0)
+ std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes, bool _isLocal)
+ : sBox(sizeBox), variableShapeParam(variableEps), radialFunctionIndex (rbfFun), XYZkdtree(0), _inUV(0), isLocal(_isLocal)
{
allCenters.resize(allNodes.size(),3);
- double tol = 6.e-2*sBox;
+ double tol = 6.e-1*sBox;
+ if (isLocal) tol = 1.e-15;
//remove duplicate vertices
//add bc nodes
@@ -111,9 +112,11 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
centers(index,1) = v1->y()/sBox;
centers(index,2) = v1->z()/sBox;
std::map<MVertex*, SVector3>::iterator itn = _normals.find(v1);
- normals(index,0) = itn->second.x();
- normals(index,1) = itn->second.y();
- normals(index,2) = itn->second.z();
+ if (itn != _normals.end()){
+ normals(index,0) = itn->second.x();
+ normals(index,1) = itn->second.y();
+ normals(index,2) = itn->second.z();
+ }
_mapV.insert(std::make_pair(v1, index));
for(std::set<MVertex *>::iterator itv2 = myNodes.begin(); itv2 !=myNodes.end() ; itv2++){
MVertex *v2 = *itv2;
@@ -125,13 +128,7 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
}
delta = 0.33*dist_min;//0.33
- radius= nbNodes/10.;
-
- matAInv.resize(nbNodes, nbNodes);
- matAInv = generateRbfMat(0,centers,centers);
- matAInv.invertInPlace();
- isLocal = false;
- extendedX.resize(3*nbNodes,3);
+ radius= 1.0/6.0; //size 1 is non dim size
Msg::Info("*****************************************");
Msg::Info("*** RBF nbNodes=%d (all=%d) ", myNodes.size(), allNodes.size());
@@ -139,7 +136,14 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
Msg::Info("*****************************************");
printNodes(myNodes);
- //exit(1);
+
+ if (!isLocal){
+ matAInv.resize(nbNodes, nbNodes);
+ matAInv = generateRbfMat(0,centers,centers);
+ matAInv.invertInPlace();
+ }
+
+ extendedX.resize(3*nbNodes,3);
}
@@ -166,6 +170,7 @@ void GRbf::buildXYZkdtree(){
void GRbf::buildOctree(double radius){
+ printf("building octree radius = %g \n", radius);
SBoundingBox3d bb;
for (int i= 0; i< nbNodes; i++)
bb += SPoint3(centers(i,0),centers(i,1), centers(i,2));
@@ -197,6 +202,7 @@ void GRbf::buildOctree(double radius){
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());
}
Octree_Delete(oct);
@@ -205,46 +211,84 @@ void GRbf::buildOctree(double radius){
}
//compute curvature from level set
-void GRbf::computeCurvature(std::map<MVertex*, SPoint3> &rbf_curv){
+
+void GRbf::curvatureRBF(const fullMatrix<double> &cntrs,
+ fullMatrix<double> &curvature){
fullMatrix<double> extX, surf, sx,sy,sz, sxx,syy,szz, sxy,sxz,syz,sLap;
+ setup_level_set(cntrs,normals,extX, surf);
- isLocal = true;
- fullMatrix<double> curvature(centers.size1(), 1);
-
- setup_level_set(centers,normals,extX, surf);
-
- evalRbfDer(1,extX,centers,surf,sx);
- evalRbfDer(2,extX,centers,surf,sy);
- evalRbfDer(3,extX,centers,surf,sz);
- evalRbfDer(11,extX,centers,surf,sxx);
- evalRbfDer(12,extX,centers,surf,sxy);
- evalRbfDer(13,extX,centers,surf,sxz);
- evalRbfDer(22,extX,centers,surf,syy);
- evalRbfDer(23,extX,centers,surf,syz);
- evalRbfDer(33,extX,centers,surf,szz);
- evalRbfDer(222,extX,centers,surf,sLap);
-
- for (int i = 0; i < centers.size1(); i++) {
+ 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;
}
+}
+
+void GRbf::computeCurvature(const fullMatrix<double> &cntrs,
+ std::map<MVertex*, double> &rbf_curv){
+
+ fullMatrix<double> curvature(cntrs.size1(), 1);
+ curvatureRBF(cntrs, curvature);
//interpolate
fullMatrix<double> curvatureAll(allCenters.size1(), 1);
- evalRbfDer(0,centers,allCenters,curvature,curvatureAll);
+ evalRbfDer(0,cntrs,allCenters,curvature,curvatureAll);
//fill rbf_curv
std::map<MVertex*, int>::iterator itm = _mapAllV.begin();
for (; itm != _mapAllV.end(); itm++){
int index = itm->second;
- SPoint3 coords(curvatureAll(index,0), 0.0, 0.0);
- rbf_curv.insert(std::make_pair(itm->first,coords));
+ rbf_curv.insert(std::make_pair(itm->first,curvatureAll(index,0)));
}
}
+void GRbf::computeLocalCurvature(const fullMatrix<double> &cntrs,
+ std::map<MVertex*, double> &rbf_curv) {
+
+ int numNodes = cntrs.size1();
+
+ if(nodesInSphere.size() == 0) buildOctree(radius);
+ fullMatrix<double> curvature(cntrs.size1(), 1);
+
+ for (int i=0;i<numNodes ; ++i){
+ std::vector<int> &pts = nodesInSphere[i];
+
+ fullMatrix<double> nodes_in_sph(pts.size(),3);
+ fullMatrix<double> LocalOper;
+
+ for (int k=0; k< pts.size(); k++){
+ nodes_in_sph(k,0) = cntrs(pts[k],0);
+ nodes_in_sph(k,1) = cntrs(pts[k],1);
+ nodes_in_sph(k,2) = cntrs(pts[k],2);
+ }
+
+ 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);
+ }
+
+ std::map<MVertex*, int>::iterator itm = _mapAllV.begin();
+ for (; itm != _mapAllV.end(); itm++){
+ int index = itm->second;
+ rbf_curv.insert(std::make_pair(itm->first,curvature(index,0)));
+ }
+
+}
+
double GRbf::evalRadialFnDer (int p, double dx, double dy, double dz, double ep){
double r2 = dx*dx+dy*dy+dz*dz; //r^2
diff --git a/Geo/GRbf.h b/Geo/GRbf.h
index e849b725d8..c3a565e77f 100644
--- a/Geo/GRbf.h
+++ b/Geo/GRbf.h
@@ -65,7 +65,7 @@ class GRbf {
GRbf (double sizeBox, int variableEps, int rbfFun,
std::map<MVertex*, SVector3> normals,
- std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes);
+ std::set<MVertex *> allNodes, std::vector<MVertex*> bcNodes, bool isLocal = false);
~GRbf();
//build octree
@@ -135,7 +135,15 @@ class GRbf {
fullMatrix<double> &D);
// Calculates the curvature of a surface at centers
- void computeCurvature(std::map<MVertex*, SPoint3> &rbf_curv);
+
+ void curvatureRBF(const fullMatrix<double> &cntrs,
+ fullMatrix<double> &curvature);
+
+ void computeCurvature(const fullMatrix<double> &cntrs,
+ std::map<MVertex*, double>&rbf_curv);
+
+ void computeLocalCurvature(const fullMatrix<double> &cntrs,
+ std::map<MVertex*, double>&rbf_curv);
//Finds the U,V,S (in the 0-level set) that are the 'num_neighbours' closest to u_eval and v_eval.
//Thus in total, we're working with '3*num_neighbours' nodes
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 38eff4782a..c46aed378d 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -503,15 +503,16 @@ double discreteEdge::curvature(double par) const{
Curvature& curvature = Curvature::getInstance();
if( !Curvature::valueAlreadyComputed() ) {
- std::cout << "Need to compute discrete curvature" << std::endl;
+ 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");
+ //curvature.writeDirectionsToPosFile("curvature_directions.pos");
+ //curvature.writeToVtkFile("curvature.vtk");
std::cout << " ... computing curvature finished" << std::endl;
}
diff --git a/benchmarks/curvature/test.geo b/benchmarks/curvature/test.geo
index 2099a02566..6135e8eab3 100644
--- a/benchmarks/curvature/test.geo
+++ b/benchmarks/curvature/test.geo
@@ -5,8 +5,7 @@
//Mesh.CharacteristicLengthFactor= 0.1;
Mesh.CharacteristicLengthFromPoints = 0;
- Mesh.CharacteristicLengthFromPoints = 0;
- Mesh.CharacteristicLengthMin = 0.01;
+ Mesh.CharacteristicLengthMin = 0.001;
Mesh.CharacteristicLengthMax = 10.00;
Mesh.LcIntegrationPrecision=1.e-5;
Mesh.MinimumCirclePoints=20;
diff --git a/benchmarks/stl/reparam.geo b/benchmarks/stl/reparam.geo
index 72d11c19d9..200571d5b4 100644
--- a/benchmarks/stl/reparam.geo
+++ b/benchmarks/stl/reparam.geo
@@ -1,8 +1,10 @@
+Mesh.RemeshParametrization=0; //(0) harmonic (1) conformal
+Mesh.RemeshAlgorithm=1; //(0) nosplit (1) automatic (2) split metis
+//Mesh.Algorithm = 6; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg)
Mesh.CharacteristicLengthFactor=0.5;
-//Geometry.Lines = 0;
-Merge "reparam_input.msh";
+Merge "reparamINPUT.msh";
CreateTopology;
Compound Line(10)={1};
diff --git a/benchmarks/stl/reparam_input.msh b/benchmarks/stl/reparamINPUT.msh
similarity index 100%
rename from benchmarks/stl/reparam_input.msh
rename to benchmarks/stl/reparamINPUT.msh
--
GitLab