From 060857d32dbc09a66f1d169df4fd035645fb34e3 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Mon, 27 Feb 2012 14:41:43 +0000
Subject: [PATCH] reintroduced harmonic mapping on square ...
---
Geo/GFaceCompound.cpp | 121 ++++++++++++++++++++++++++++++------------
Geo/GFaceCompound.h | 22 ++++++--
Geo/GModelIO_Geo.cpp | 16 +++---
Mesh/meshGFace.cpp | 3 +-
4 files changed, 115 insertions(+), 47 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 0f525255e4..93ac515eb7 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -683,7 +683,6 @@ 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; // 1 MQ RBF , 0 GA
double sizeBox = getSizeH();
@@ -695,17 +694,8 @@ bool GFaceCompound::parametrize() const
//_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);
-
- _rbf->solveHarmonicMap(Oper, _ordered, _coords, coordinates);
-
- //_rbf->computeCurvature(coordinates);
- //printStuff();
- //exit(1);
- double t2 = Cpu();
- printf("param RBF in %g sec \n", t2-t1);
+ _rbf->solveHarmonicMap(Oper, _ordered, _coords, coordinates);
}
buildOct();
@@ -750,6 +740,23 @@ void GFaceCompound::getBoundingEdges()
l_edges.push_back(*it);
(*it)->addFace(this);
}
+ if (_V0.size() && _U1.size() && _V1.size()){
+ it = _V0.begin();
+ for( ; it != _V0.end() ; ++it){
+ l_edges.push_back(*it);
+ (*it)->addFace(this);
+ }
+ it = _U1.begin();
+ for( ; it != _U1.end() ; ++it){
+ l_edges.push_back(*it);
+ (*it)->addFace(this);
+ }
+ it = _V1.begin();
+ for( ; it != _V1.end() ; ++it){
+ l_edges.push_back(*it);
+ (*it)->addFace(this);
+ }
+ }
std::list<GEdge*> loop;
computeALoop(_unique, loop);
while(!_unique.empty()) computeALoop(_unique,loop);
@@ -976,6 +983,39 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
nbSplit = 0;
fillTris.clear();
}
+GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
+ std::list<GEdge*> &U0, std::list<GEdge*> &V0,
+ std::list<GEdge*> &U1, std::list<GEdge*> &V1,
+ typeOfMapping mpg,
+ int allowPartition,
+ linearSystem<double>* lsys)
+ : GFace(m, tag), _compound(compound), oct(0),
+ _U0(U0), _V0(V0), _U1(U1), _V1(V1),
+ _mapping(mpg), _allowPartition(allowPartition), _lsys(lsys)
+{
+ ONE = new simpleFunction<double>(1.0);
+ MONE = new simpleFunction<double>(-1.0);
+
+ for(std::list<GFace*>::iterator it = _compound.begin(); it != _compound.end(); ++it){
+ if(!(*it)){
+ Msg::Error("Incorrect face in compound surface %d\n", tag);
+ Msg::Exit(1);
+ }
+ }
+
+ getBoundingEdges();
+ _type = UNITCIRCLE;
+ if (_mapping == HARMONICPLANE){
+ _mapping = HARMONIC;
+ _type = MEANPLANE;
+ }
+ else if(_mapping == HARMONIC && _U0.size() && _V0.size() && _U1.size() && _V1.size()) {
+ _type = SQUARE;
+ }
+
+ nbSplit = 0;
+ fillTris.clear();
+}
GFaceCompound::~GFaceCompound()
{
@@ -1101,39 +1141,50 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom) const
else if(step == ITERV) myAssembler.fixVertex(v, 0, 1, sin(theta));
}
}
+ else if(_type == SQUARE){
+ if(step == ITERU){
+ std::list<GEdge*>::const_iterator it = _U0.begin();
+ for( ; it != _U0.end(); ++it){
+ fixEdgeToValue(*it, 0.0, myAssembler);
+ }
+ it = _U1.begin();
+ for( ; it != _U1.end(); ++it){
+ fixEdgeToValue(*it, 1.0, myAssembler);
+ }
+ }
+ else if(step == ITERV){
+ std::list<GEdge*>::const_iterator it = _V0.begin();
+ for( ; it != _V0.end(); ++it){
+ fixEdgeToValue(*it, 0.0, myAssembler);
+ }
+ it = _V1.begin();
+ for( ; it != _V1.end(); ++it){
+ fixEdgeToValue(*it, 1.0, myAssembler);
+ }
+ }
+ }
else if (_type == MEANPLANE){
-
- FILE * f3 = fopen("PTS.pos","w");
- fprintf(f3, "View \"\"{\n");
- std::vector<SPoint3> points, pointsProj, pointsUV;
- SPoint3 ptCG(0.0, 0.0, 0.0);
- for(unsigned int i = 0; i < _ordered.size(); i++){
+ std::vector<SPoint3> points, pointsProj, pointsUV;
+ SPoint3 ptCG(0.0, 0.0, 0.0);
+ for(unsigned int i = 0; i < _ordered.size(); i++){
MVertex *v = _ordered[i];
points.push_back(SPoint3(v->x(), v->y(), v->z()));
- fprintf(f3, "SP(%g,%g,%g){%g};\n",
- v->x(), v->y(), v->z(),
- (double)v->getNum());
ptCG[0] += v->x();
ptCG[1] += v->y();
ptCG[2] += v->z();
}
- fprintf(f3,"};\n");
- fclose(f3);
-
- ptCG /= points.size();
- mean_plane meanPlane;
- computeMeanPlaneSimple(points, meanPlane);
- projectPointsToPlane(points, pointsProj, meanPlane);
- transformPointsIntoOrthoBasis(pointsProj, pointsUV, ptCG, meanPlane);
-
- for(unsigned int i = 0; i < pointsUV.size(); i++){
- MVertex *v = _ordered[i];
- if(step == ITERU) myAssembler.fixVertex(v, 0, 1, pointsUV[i][0]);
- else if(step == ITERV) myAssembler.fixVertex(v, 0, 1, pointsUV[i][1]);
- }
-
+ ptCG /= points.size();
+ mean_plane meanPlane;
+ computeMeanPlaneSimple(points, meanPlane);
+ projectPointsToPlane(points, pointsProj, meanPlane);
+ transformPointsIntoOrthoBasis(pointsProj, pointsUV, ptCG, meanPlane);
+ for(unsigned int i = 0; i < pointsUV.size(); i++){
+ MVertex *v = _ordered[i];
+ if(step == ITERU) myAssembler.fixVertex(v, 0, 1, pointsUV[i][0]);
+ else if(step == ITERV) myAssembler.fixVertex(v, 0, 1, pointsUV[i][1]);
+ }
}
else{
Msg::Error("Unknown type of parametrization");
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index d666100b88..f3d898e4c3 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -57,7 +57,7 @@ class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
typedef enum {HARMONIC=1,CONFORMAL=2, RBF=3, HARMONICPLANE=4, CONVEXCOMBINATION=5} typeOfMapping;
- typedef enum {UNITCIRCLE, MEANPLANE} typeOfIsomorphism;
+ typedef enum {UNITCIRCLE, MEANPLANE, SQUARE} typeOfIsomorphism;
void computeNormals(std::map<MVertex*, SVector3> &normals) const;
protected:
mutable std::set<MVertex *> ov;
@@ -65,7 +65,7 @@ class GFaceCompound : public GFace {
simpleFunction<double> *ONE;
simpleFunction<double> *MONE;
std::list<GFace*> _compound;
- std::list<GEdge*> _U0;
+ std::list<GEdge*> _U0, _V0, _U1, _V1;
std::list<std::list<GEdge*> > _interior_loops;
mutable int nbT;
mutable GFaceCompoundTriangle *_gfct;
@@ -113,6 +113,12 @@ class GFaceCompound : public GFace {
std::list<GEdge*> &U0, typeOfMapping typ = HARMONIC,
int allowPartition=1,
linearSystem<double>* lsys =0);
+ GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
+ std::list<GEdge*> &U0, std::list<GEdge*> &V0,
+ std::list<GEdge*> &U1, std::list<GEdge*> &V1,
+ typeOfMapping typ = HARMONIC,
+ int allowPartition=1,
+ linearSystem<double>* lsys =0);
virtual ~GFaceCompound();
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
@@ -151,7 +157,7 @@ class GFaceCompound : public GFace {
public:
typedef enum {ITERU=0,ITERV=1,ITERD=2} iterationStep;
typedef enum {HARMONIC=1,CONFORMAL=2, RBF=3, HARMONICPLANE=4, CONVEXCOMBINATION=5} typeOfMapping;
- typedef enum {UNITCIRCLE, MEANPLANE} typeOfIsomorphism;
+ typedef enum {UNITCIRCLE, MEANPLANE, SQUARE} typeOfIsomorphism;
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfMapping typ = HARMONIC,
int allowPartition=1,
@@ -160,6 +166,16 @@ class GFaceCompound : public GFace {
{
Msg::Error("Gmsh has to be compiled with solver support to use GFaceCompounds");
}
+ GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
+ std::list<GEdge*> &U0, std::list<GEdge*> &V0,
+ std::list<GEdge*> &U1, std::list<GEdge*> &V1,
+ typeOfMapping typ = HARMONIC,
+ int allowPartition=1,
+ linearSystem<double>* lsys =0)
+ : GFace(m, tag)
+ {
+ Msg::Error("Gmsh has to be compiled with solver support to use GFaceCompounds");
+ }
virtual ~GFaceCompound() {}
GPoint point(double par1, double par2) const { return GPoint(); }
Pair<SVector3, SVector3> firstDer(const SPoint2 ¶m) const
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 985df492a5..a0f921a258 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -182,22 +182,22 @@ int GModel::importGEOInternals()
if(gf)
comp.push_back(gf);
}
- std::list<GEdge*> U0;
+ std::list<GEdge*> b[4];
for(int j = 0; j < 4; j++){
- for(unsigned int k = 0; k < s->compoundBoundary[j].size(); k++){
- GEdge *ge = getEdgeByTag(s->compoundBoundary[j][k]);
- if(ge) U0.push_back(ge);
- }
- }
+ for(unsigned int k = 0; k < s->compoundBoundary[j].size(); k++){
+ GEdge *ge = getEdgeByTag(s->compoundBoundary[j][k]);
+ if(ge) b[j].push_back(ge);
+ }
+ }
int param = CTX::instance()->mesh.remeshParam;
GFaceCompound::typeOfMapping typ = GFaceCompound::HARMONIC;
if (param == 1) typ = GFaceCompound::CONFORMAL;
if (param == 2) typ = GFaceCompound::RBF;
if (param == 3) typ = GFaceCompound::HARMONICPLANE;
- if (param == 4) typ = GFaceCompound::CONVEXCOMBINATION;
+ if (param == 4) typ = GFaceCompound::CONVEXCOMBINATION;
int algo = CTX::instance()->mesh.remeshAlgo;
- f = new GFaceCompound(this, s->Num, comp, U0, typ, algo);
+ f = new GFaceCompound(this, s->Num, comp, b[0], b[1], b[2], b[3], typ, algo);
f->meshAttributes.recombine = s->Recombine;
f->meshAttributes.recombineAngle = s->RecombineAngle;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e36d70aa73..6b7d64a66d 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1749,10 +1749,11 @@ void deMeshGFace::operator() (GFace *gf)
}
// for debugging, change value from -1 to -100;
-int debugSurface = -1; //-1;
+int debugSurface = -100; //-1;
void meshGFace::operator() (GFace *gf)
{
+
gf->model()->setCurrentMeshEntity(gf);
if(debugSurface >= 0 && gf->tag() != debugSurface){
--
GitLab