From e6496ee8cc25d85d9c41d7247aa06d80e0942d86 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Thu, 18 Aug 2011 15:11:58 +0000
Subject: [PATCH] Explicit not yet implemented for incomp (should add
artificial compressibility)
---
Geo/GFaceCompound.cpp | 26 ++++++++++++++++++++++----
Geo/GRbf.cpp | 13 +++++++------
Solver/function.cpp | 2 +-
3 files changed, 30 insertions(+), 11 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index b6f7b69c36..8922a80eda 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -83,6 +83,19 @@ int intersection_segments (SPoint3 &p1, SPoint3 &p2,
x[1] >= 0.0 && x[1] <= 1.);
}
}
+static void printBound(std::vector<MVertex*> &l, int tag){
+ //print boundary
+ char name[256];
+ sprintf(name, "myBOUND%d.pos", tag);
+ FILE * xyz = fopen(name,"w");
+ fprintf(xyz,"View \"\"{\n");
+ for(int i = 0; i < l.size(); i++){
+ MVertex *v = l[i];
+ fprintf(xyz,"SP(%g,%g,%g){%d};\n", v->x(), v->y(), v->z(), i);
+ }
+ fprintf(xyz,"};\n");
+ fclose(xyz);
+}
static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
std::vector<double> &coord)
@@ -144,6 +157,7 @@ static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
}
if(!found) return false;
}
+
return true;
}
@@ -628,7 +642,8 @@ bool GFaceCompound::parametrize() const
if(allNodes.empty()) buildAllNodes();
bool success = orderVertices(_U0, _ordered, _coords);
- if(!success) Msg::Error("Could not order vertices on boundary");
+ printBound(_ordered, this->tag());
+ if(!success) {Msg::Error("Could not order vertices on boundary");exit(1);}
// Laplace parametrization
if (_mapping == HARMONIC){
@@ -636,6 +651,8 @@ bool GFaceCompound::parametrize() const
fillNeumannBCS();
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
+ //parametrize(ITERU,CONVEXCOMBINATION);
+ //parametrize(ITERV,CONVEXCOMBINATION);
}
// Multiscale Laplace parametrization
else if (_mapping == MULTISCALE){
@@ -663,9 +680,9 @@ bool GFaceCompound::parametrize() const
// Radial-Basis Function parametrization
else if (_mapping == RBF){
Msg::Debug("Parametrizing surface %d with 'RBF' ", tag());
-
+ double t1 = Cpu();
int variableEps = 0;
- int radFunInd = 1; //MQ RBF
+ int radFunInd = 1; // 1 MQ RBF , 0 GA
double sizeBox = getSizeH();
fullMatrix<double> Oper(3*allNodes.size(),3*allNodes.size());
@@ -683,7 +700,8 @@ bool GFaceCompound::parametrize() const
//_rbf->computeCurvature(coordinates);
//printStuff();
//exit(1);
-
+ double t2 = Cpu();
+ printf("param RBF in %g sec \n", t2-t1);
}
buildOct();
diff --git a/Geo/GRbf.cpp b/Geo/GRbf.cpp
index 94195882f6..52f9610d75 100644
--- a/Geo/GRbf.cpp
+++ b/Geo/GRbf.cpp
@@ -74,13 +74,13 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
{
allCenters.resize(allNodes.size(),3);
- double tol = 4.e-2*sBox;
+ double tol = 6.e-2*sBox;
//remove duplicate vertices
//add bc nodes
for(unsigned int i = 0; i < bcNodes.size(); i++){
myNodes.insert(bcNodes[i]);
- //if (bcNodes.size() > 20) i+=3;
+ //if (bcNodes.size() > 20) i+=2;
}
//then create Mvertex position
@@ -104,6 +104,7 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
normals.resize(nbNodes,3);
int index = 0;
double dist_min = 1.e6;
+ double dist_max = 1.e-6;
for(std::set<MVertex *>::iterator itv = myNodes.begin(); itv !=myNodes.end() ; ++itv){
MVertex *v1 = *itv;
centers(index,0) = v1->x()/sBox;
@@ -118,12 +119,12 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
MVertex *v2 = *itv2;
double dist = sqrt((v1->x()-v2->x())*(v1->x()-v2->x())+(v1->y()-v2->y())*(v1->y()-v2->y())
+(v1->z()-v2->z())*(v1->z()-v2->z()))/sBox;
- if (dist<dist_min && *itv != *itv2 && dist > 1.e-6) dist_min = dist;
+ if (dist<dist_min && *itv != *itv2 && dist > 1.e-5) dist_min = dist;
}
index++;
}
- delta = 0.33*dist_min;//0.6
+ delta = 0.33*dist_min;//0.33
radius= nbNodes/10.;
matAInv.resize(nbNodes, nbNodes);
@@ -134,7 +135,7 @@ GRbf::GRbf (double sizeBox, int variableEps, int rbfFun, std::map<MVertex*, SVec
Msg::Info("*****************************************");
Msg::Info("*** RBF nbNodes=%d (all=%d) ", myNodes.size(), allNodes.size());
- Msg::Info("*** RBF rad=%g dist_min =%g ", radius, dist_min);
+ Msg::Info("*** RBF rad=%g dist_min =%g", radius, dist_min);
Msg::Info("*****************************************");
printNodes(myNodes);
@@ -843,7 +844,7 @@ bool GRbf::UVStoXYZ(const double u_eval, const double v_eval,
for (int j = i+1; j < num_neighbours; j++){
double dist = sqrt((UV(index[i],0)-UV(index[j],0))*(UV(index[i],0)-UV(index[j],0))+
(UV(index[i],1)-UV(index[j],1))*(UV(index[i],1)-UV(index[j],1)));
- if (dist<dist_min && dist > 1.e-6) dist_min = dist;
+ if (dist<dist_min && dist > 1.e-5) dist_min = dist;
}
}
delete [] index;
diff --git a/Solver/function.cpp b/Solver/function.cpp
index 5f9f3d9a71..f589a260a7 100644
--- a/Solver/function.cpp
+++ b/Solver/function.cpp
@@ -467,7 +467,6 @@ class functionLevelsetSmooth : public function {
double _valMin, _valPlus, _E;
void call(dataCacheMap *m, fullMatrix<double> &val)
{
-
double ivalPlus = 1./_valPlus;
double ivalMin = 1./_valMin;
for (int i = 0; i < val.size1(); i++)
@@ -491,6 +490,7 @@ class functionLevelsetSmooth : public function {
}
functionLevelsetSmooth(const function *f0, const double valMin, const double valPlus, const double E) : function(f0->getNbCol())
{
+ printf("Levelset bandwidth is E = %g \n", E);
setArgument (_f0, f0);
_valMin = valMin;
_valPlus = valPlus;
--
GitLab