From c73a077d697465c0b4eeeb475c6693f108debca2 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 12 Aug 2009 09:24:14 +0000
Subject: [PATCH] *** empty log message ***
---
Geo/GEdgeCompound.cpp | 23 +++++++++
Geo/GEdgeCompound.h | 4 +-
Mesh/meshGFaceBDS.cpp | 88 ++++++++++++++++++++++++++++++---
Numeric/gmshLinearSystemCSR.cpp | 2 +-
4 files changed, 108 insertions(+), 9 deletions(-)
diff --git a/Geo/GEdgeCompound.cpp b/Geo/GEdgeCompound.cpp
index 1e40b3984d..2a4c65742d 100644
--- a/Geo/GEdgeCompound.cpp
+++ b/Geo/GEdgeCompound.cpp
@@ -195,6 +195,29 @@ void GEdgeCompound::getLocalParameter ( const double &t,
}
}
+// give a GEdge and
+void GEdgeCompound::getCompoundParameter ( GEdge *ge,
+ const double &tLoc,
+ double &t) const
+{
+
+ for (int iEdge = 0; iEdge < (int)_compound.size(); iEdge++){
+ //printf("iEdge=%d tmin=%g\n",iEdge,_pars[iEdge]);
+ if (ge == _compound[iEdge]){
+ double tmin = _pars[iEdge];
+ double tmax = _pars[iEdge+1];
+ Range<double> b = _compound[iEdge]->parBounds(0);
+ t = _orientation[iEdge] ?
+ tmin + (tLoc - b.low())/(b.high()-b.low()) * (tmax-tmin):
+ tmax - (tLoc - b.low())/(b.high()-b.low()) * (tmax-tmin);
+ //printf("bhigh=%g, blow=%g, global t=%g , tLoc=%g ,iEdge=%d\n",b.high(), b.low(), t,tLoc,iEdge);
+ return;
+ }
+ }
+}
+
+
+
void GEdgeCompound::parametrize()
{
_pars.push_back(0.0);
diff --git a/Geo/GEdgeCompound.h b/Geo/GEdgeCompound.h
index 807b717528..0e5a272b6e 100644
--- a/Geo/GEdgeCompound.h
+++ b/Geo/GEdgeCompound.h
@@ -24,7 +24,9 @@ public:
void getLocalParameter ( const double &t,
int &iEdge,
double & tLoc) const;
-
+ void getCompoundParameter ( GEdge *ge,
+ const double &tLoc,
+ double &t) const;
GEdgeCompound(GModel *m, int tag, std::vector<GEdge*> &compound);
virtual ~GEdgeCompound();
Range<double> parBounds(int i) const;
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index e2b293f755..92c877b791 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -33,10 +33,10 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
}
inline double computeEdgeLinearLength(BDS_Point *p1,
- BDS_Point *p2,
- GFace *f,
- double SCALINGU,
- double SCALINGV)
+ BDS_Point *p2,
+ GFace *f,
+ double SCALINGU,
+ double SCALINGV)
{
GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
@@ -56,6 +56,78 @@ inline double computeEdgeLinearLength(BDS_Point *p1,
return l1 + l2;
}
+inline double computeEdgeLinearLength_new(BDS_Point *p1,
+ BDS_Point *p2,
+ GFace *f,
+ double SCALINGU,
+ double SCALINGV)
+{
+ const int nbSb = 10;
+ GPoint GP[nbSb-1];
+
+ for (int i=1;i<nbSb;i++){
+ double xi = (double)i/nbSb;
+ GP[i-1] = f->point(SPoint2(((1-xi) * p1->u + xi * p2->u) * SCALINGU,
+ ((1-xi) * p1->v + xi * p2->v) * SCALINGV));
+ if (!GP[i-1].succeeded())
+ return computeEdgeLinearLength(p1,p2);
+ }
+ double l = 0;
+ for (int i=0;i<nbSb;i++){
+ if (i == 0){
+ const double dx1 = p1->X - GP[0].x();
+ const double dy1 = p1->Y - GP[0].y();
+ const double dz1 = p1->Z - GP[0].z();
+ l+= sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
+ }
+ else if (i == nbSb-1){
+ const double dx1 = p2->X - GP[nbSb-1].x();
+ const double dy1 = p2->Y - GP[nbSb-1].y();
+ const double dz1 = p2->Z - GP[nbSb-1].z();
+ l+= sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
+ }
+ else{
+ const double dx1 = GP[i].x() - GP[i-1].x();
+ const double dy1 = GP[i].y() - GP[i-1].y();
+ const double dz1 = GP[i].z() - GP[i-1].z();
+ l+= sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
+ }
+ }
+ return l;
+}
+
+
+
+inline double computeEdgeMiddleCoord_new(BDS_Point *p1,
+ BDS_Point *p2,
+ GFace *f,
+ double SCALINGU,
+ double SCALINGV)
+{
+ const int nbSb = 3;
+ double L = computeEdgeLinearLength(p1,p2);
+ GPoint GP[nbSb];
+ for (int i=1;i<nbSb;i++){
+ double xi = (double)i/nbSb;
+ GP[i-1] = f->point(SPoint2(((1-xi) * p1->u + xi * p2->u) * SCALINGU,
+ ((1-xi) * p1->v + xi * p2->v) * SCALINGV));
+ if (!GP[i-1].succeeded())
+ return 0.5;
+ const double dx1 = p1->X - GP[i-1].x();
+ const double dy1 = p1->Y - GP[i-1].y();
+ const double dz1 = p1->Z - GP[i-1].z();
+ double LPLUS = sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
+ if (LPLUS > L*.5){
+ double XIMINUS, XIPLUS,LPLUS,LMINUS;
+ if (i==1){
+ XIMINUS=0;
+ }
+ return XIMINUS + (LPLUS - L*.5)/(LPLUS-LMINUS)/(nbSb-1);
+ }
+ }
+ return 0.5;
+}
+
inline double computeEdgeMiddleCoord(BDS_Point *p1,
BDS_Point *p2,
GFace *f,
@@ -431,9 +503,9 @@ void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
if (!(*it)->deleted){
double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
if ((*it)->numfaces() == 2 && (lone > MAXE_)){
- //const double coord = 0.5;
- const double coord = computeEdgeMiddleCoord((*it)->p1, (*it)->p2, gf,
- m.scalingU, m.scalingV);
+ const double coord = 0.5;
+ //const double coord = computeEdgeMiddleCoord((*it)->p1, (*it)->p2, gf,
+ // m.scalingU, m.scalingV);
BDS_Point *mid;
GPoint gpp = gf->point(m.scalingU*(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u),
@@ -478,6 +550,8 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
BDS_Edge *e = edges[i].second;
if (!e->deleted){
const double coord = 0.5;
+ // const double coord = computeEdgeMiddleCoord((*it)->p1, (*it)->p2, gf,
+ // m.scalingU, m.scalingV);
BDS_Point *mid ;
double U = coord * e->p1->u + (1 - coord) * e->p2->u;
double V = coord * e->p1->v + (1 - coord) * e->p2->v;
diff --git a/Numeric/gmshLinearSystemCSR.cpp b/Numeric/gmshLinearSystemCSR.cpp
index d80f80f8ca..01f9c03b7b 100644
--- a/Numeric/gmshLinearSystemCSR.cpp
+++ b/Numeric/gmshLinearSystemCSR.cpp
@@ -291,7 +291,7 @@ void sortColumns (int NbLines,
template<>
int gmshLinearSystemCSRGmm<double> :: systemSolve()
{
- if(!sorted)
+ if (!sorted)
sortColumns(_b->size(),
CSRList_Nbr(a_),
(INDEX_TYPE *) ptr_->array,
--
GitLab