*** empty log message ***

parent 1c7caca8
......@@ -1198,7 +1198,6 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
GPoint gp = gf->point(U * scalingU, V * scalingV);
if (!gp.succeeded()){
// printf ("iha\n");
return false;
}
const double oldX = p->X;
......
......@@ -77,16 +77,16 @@ static double LC_MVertex_CURV(GEntity *ge, double U, double V)
double Crv = 0;
switch(ge->dim()){
case 0:
Crv = max_edge_curvature((const GVertex *)ge);
Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
//Crv = max_surf_curvature((const GVertex *)ge);
// Crv = max_edge_curvature((const GVertex *)ge);
// Crv = std::max(max_surf_curvature((const GVertex *)ge), Crv);
Crv = max_surf_curvature((const GVertex *)ge);
break;
case 1:
{
GEdge *ged = (GEdge *)ge;
Crv = ged->curvature(U);
Crv = std::max(Crv, max_surf_curvature(ged, U));
//Crv = max_surf_curvature(ged, U);
// Crv = ged->curvature(U)*2;
// Crv = std::max(Crv, max_surf_curvature(ged, U));
Crv = max_surf_curvature(ged, U);
}
break;
case 2:
......
......@@ -496,14 +496,15 @@ static int ConvertDListToTriangles(DocRecord *doc)
{
// on suppose que n >= 3. gPointArray est suppose OK.
STriangle *striangle;
// STriangle *striangle;
int n, i, j;
int count = 0, count2 = 0;
PointNumero aa, bb, cc;
PointRecord *gPointArray = doc->points;
n = doc->numPoints;
striangle = (STriangle *) Malloc(n * sizeof(STriangle));
STriangle *striangle = (STriangle *) Malloc(n * sizeof(STriangle));
count2 = CountPointsOnHull(n, doc->points);
// nombre de triangles que l'on doit obtenir
......@@ -516,8 +517,6 @@ static int ConvertDListToTriangles(DocRecord *doc)
// de points (t_length)
striangle[i].t = ConvertDlistToArray(&gPointArray[i].adjacent,
&striangle[i].t_length);
striangle[i].info = NULL;
striangle[i].info_length = 0;
}
// on balaye les noeuds de gauche a droite -> on cree les triangles
......@@ -536,11 +535,8 @@ static int ConvertDListToTriangles(DocRecord *doc)
}
}
}
for(i = 0; i < n; i++)
Free(striangle[i].t);
Free(striangle);
doc->numTriangles = count2;
doc->striangle = striangle;
return 1;
}
......@@ -563,7 +559,7 @@ static void RemoveAllDList(int n, PointRecord *pPointArray)
}
DocRecord::DocRecord(int n)
: numPoints(n), points(NULL), numTriangles(0), triangles(NULL)
: numPoints(n), points(NULL), striangle(0),numTriangles(0), triangles(NULL)
{
if(numPoints)
points = (PointRecord*)Malloc(numPoints * sizeof(PointRecord));
......@@ -573,6 +569,11 @@ DocRecord::~DocRecord()
{
Free(points);
Free(triangles);
if (striangle){
for(int i = 0; i < numPoints; i++)
Free(striangle[i].t);
Free(striangle);
}
}
void DocRecord::MakeMeshWithPoints()
......@@ -581,4 +582,17 @@ void DocRecord::MakeMeshWithPoints()
BuildDelaunay(this);
ConvertDListToTriangles(this);
RemoveAllDList(numPoints, points);
for(int i = 0; i < numPoints; i++)
Free(striangle[i].t);
Free(striangle);
striangle = 0;
}
void DocRecord::MakeVoronoi(){
if(numPoints < 3) return;
BuildDelaunay(this);
// TagPointsOnHull(this);
ConvertDListToTriangles(this);
RemoveAllDList(numPoints, points);
}
......@@ -26,14 +26,8 @@ struct _CDLIST{
};
typedef struct{
PointNumero search;
PointNumero already;
}HalfTriangle;
typedef struct{
HalfTriangle *info;
PointNumero *t;
int t_length, info_length;
int t_length;
}STriangle;
typedef struct{
......@@ -50,15 +44,22 @@ typedef struct{
PointNumero a, b, c;
}Triangle;
class DocRecord{
public:
int numPoints; // number of points
PointRecord *points; // points to triangulate
int numTriangles; // number of triangles
Triangle *triangles; // 2D results
STriangle *striangle; // connected points
DocRecord(int n);
double & x(int i)
{return points[i].where.v;}
double & y(int i)
{return points[i].where.h;}
~DocRecord();
void MakeMeshWithPoints();
void MakeVoronoi();
};
#endif
......@@ -149,6 +149,7 @@ gmshSmoothHighOrder${OBJEXT}: gmshSmoothHighOrder.cpp ../Geo/MLine.h \
../Numeric/GmshMatrix.h ../Numeric/gmshElasticity.h \
../Numeric/gmshTermOfFormulation.h ../Numeric/GmshMatrix.h \
../Numeric/gmshLinearSystemGmm.h ../Numeric/gmshLinearSystem.h \
../Numeric/gmshLinearSystemCSR.h ../Numeric/gmshLinearSystem.h \
../Common/Context.h ../Geo/CGNSOptions.h ../Mesh/meshPartitionOptions.h
meshGEdge${OBJEXT}: meshGEdge.cpp ../Geo/GModel.h ../Geo/GVertex.h \
../Geo/GEntity.h ../Geo/Range.h ../Geo/SPoint3.h \
......
......@@ -21,6 +21,7 @@
#include "gmshLaplace.h"
#include "gmshElasticity.h"
#include "gmshLinearSystemGmm.h"
#include "gmshLinearSystemCSR.h"
#include "GFace.h"
#include "GRegion.h"
#include "Context.h"
......@@ -428,9 +429,9 @@ void gmshHighOrderSmoother::computeMetricVector (GFace *gf,
void gmshHighOrderSmoother::smooth_metric ( std::vector<MElement*> & all, GFace *gf)
{
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
// lsys->setNoisy(1);
lsys->setNoisy(1);
lsys->setGmres(1);
lsys->setPrec(1.e-7);
lsys->setPrec(1.e-9);
gmshAssembler<double> myAssembler(lsys);
gmshElasticityTerm El(0, 1.0, .333, getTag());
......@@ -516,7 +517,7 @@ void gmshHighOrderSmoother::smooth_metric ( std::vector<MElement*> & all, GFace
double dx2 = smooth_metric_ ( v,gf, myAssembler, verticesToMove,El);
printf(" dx2 = %12.5E\n",dx2);
if (fabs(dx2-dx) < 1.e-4 * dx0)break;
if (iter++ > 10)break;
if (iter++ > 2)break;
dx = dx2;
}
......@@ -615,10 +616,10 @@ double gmshHighOrderSmoother::smooth_metric_ ( std::vector<MElement*> & v,
void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
{
gmshLinearSystemGmm<double> *lsys = new gmshLinearSystemGmm<double>;
gmshLinearSystemCSRGmm<double> *lsys = new gmshLinearSystemCSRGmm<double>;
lsys->setNoisy(1);
lsys->setGmres(1);
lsys->setPrec(1.e-5);
lsys->setPrec(5.e-8);
gmshAssembler<double> myAssembler(lsys);
gmshElasticityTerm El(0, 1.0, .333, getTag());
......@@ -632,7 +633,7 @@ void gmshHighOrderSmoother::smooth ( std::vector<MElement*> & all)
if (!v.size()) return;
const int nbLayers = 2;
const int nbLayers = 1;
for (int i=0;i<nbLayers;i++){
addOneLayer ( all, v, layer);
v.insert(v.end(),layer.begin(),layer.end());
......
......@@ -633,7 +633,7 @@ static bool gmsh2DMeshGenerator(GFace *gf, int RECUR_ITER,
if(RECUR_ITER > 0)
Msg::Warning(":-) Gmsh was able to recover all edges after %d iterations", RECUR_ITER);
Msg::Info("Boundary Edges recovered for surface %d", gf->tag());
Msg::Debug("Boundary Edges recovered for surface %d", gf->tag());
// Look for an edge that is on the boundary for which one of the two
// neighbors has a negative number node. The other triangle is
......@@ -1077,8 +1077,8 @@ static bool gmsh2DMeshGeneratorPeriodic(GFace *gf, bool debug = true)
// Buid a BDS_Mesh structure that is convenient for doing the actual
// meshing procedure
BDS_Mesh *m = new BDS_Mesh;
m->scalingU = fabs(du);
m->scalingV = fabs(dv);
m->scalingU = 1;//fabs(du);
m->scalingV = 1;//fabs(dv);
std::vector<std::vector<BDS_Point*> > edgeLoops_BDS;
SBoundingBox3d bbox;
int nbPointsTotal = 0;
......@@ -1449,8 +1449,9 @@ void orientMeshGFace::operator()(GFace *gf)
{
gf->model()->setCurrentMeshEntity(gf);
return;
if(gf->geomType() == GEntity::ProjectionFace) return;
if(gf->geomType() == GEntity::CompoundSurface)return;
// in old versions we did not reorient transfinite surface meshes;
// we could add the following to provide backward compatibility:
// if(gf->meshAttributes.Method == MESH_TRANSFINITE) return;
......
......@@ -100,20 +100,41 @@ double NewGetLc(BDS_Point *p)
std::min(p->lc(), p->lcBGM()) : p->lcBGM();
}
static double correctLC_ (BDS_Point *p1,BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV){
double l1 = NewGetLc(p1);
double l2 = NewGetLc(p2);
double l = 0.5*(l1+l2);
// printf(" %g %g -- %g %g\n",SCALINGU,SCALINGV,l1,l2);
if(CTX::instance()->mesh.lcFromCurvature)
{
// GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
// 0.5 * (p1->v + p2->v) * SCALINGV));
// double l3 = BGM_MeshSize(f,GP.u(),GP.v(),GP.x(),GP.y(),GP.z());
double l3=l2;
double lcmin = std::min(std::min(l1,l2),l3);
l1 = std::min(lcmin*1.2,l1);
l2 = std::min(lcmin*1.2,l2);
l3 = std::min(lcmin*1.2,l3);
l = (l1+l2+l3)/3.0;
}
return l;
}
double NewGetLc(BDS_Edge *e, GFace *f, double SCALINGU, double SCALINGV)
{
double linearLength = computeEdgeLinearLength(e, f, SCALINGU, SCALINGV);
double l1 = NewGetLc(e->p1);
double l2 = NewGetLc(e->p2);
return 2 * linearLength / (l1 + l2);
double l = correctLC_ (e->p1,e->p2,f, SCALINGU, SCALINGV);
return linearLength / l;
}
double NewGetLc(BDS_Point *p1,BDS_Point *p2, GFace *f, double su, double sv)
{
double linearLength = computeEdgeLinearLength(p1, p2, f, su, sv);
double l1 = NewGetLc(p1);
double l2 = NewGetLc(p2);
return 2 * linearLength / (l1 + l2);
double l = correctLC_ (p1,p2,f, su, sv);
return linearLength / l;
}
void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
......@@ -415,8 +436,8 @@ void splitEdgePassUnsorted(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
m.scalingU, m.scalingV);
BDS_Point *mid;
GPoint gpp = gf->point(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v);
GPoint gpp = gf->point(m.scalingU*(coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u),
m.scalingV*(coord * (*it)->p1->v + (1 - coord) * (*it)->p2->v));
if (gpp.succeeded()){
mid = m.add_point(++m.MAXPOINTNUMBER,
coord * (*it)->p1->u + (1 - coord) * (*it)->p2->u,
......@@ -461,7 +482,7 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
double U = coord * e->p1->u + (1 - coord) * e->p2->u;
double V = coord * e->p1->v + (1 - coord) * e->p2->v;
GPoint gpp = gf->point(U,V);
GPoint gpp = gf->point(m.scalingU*U,m.scalingV*V);
if (gpp.succeeded()){
mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(),gpp.y(),gpp.z());
mid->u = U;
......
......@@ -69,7 +69,9 @@ void circumCenterMetric(double *pa, double *pb, double *pc,
d * (pa[1] * pa[1] - pc[1] * pc[1]) +
2. * b * (pa[0] * pa[1] - pc[0] * pc[1]);
sys2x2(sys, rhs, x);
if (!sys2x2(sys, rhs, x)){
// printf("%g %g %g\n",a,b,d);
}
Radius2 =
(x[0] - pa[0]) * (x[0] - pa[0]) * a +
......@@ -618,15 +620,15 @@ static void insertAPoint(GFace *gf,
uv[0] * vSizes [ptin->tri()->getVertex(1)->getNum()] +
uv[1] * vSizes [ptin->tri()->getVertex(2)->getNum()]);
double lc = BGM_MeshSize(gf,center[0],center[1],p.x(),p.y(),p.z());
SMetric3 metr = BGM_MeshMetric(gf,center[0],center[1],p.x(),p.y(),p.z());
//SMetric3 metr = BGM_MeshMetric(gf,center[0],center[1],p.x(),p.y(),p.z());
// vMetricsBGM.push_back(metr);
vSizesBGM.push_back(lc);
vSizes.push_back(lc1);
Us.push_back(center[0]);
Vs.push_back(center[1]);
if (!insertVertex(gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric) || !p.succeeded()) {
if (!p.succeeded() || !insertVertex(gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric) ) {
// Msg::Debug("2D Delaunay : a cavity is not star shaped");
AllTris.erase(it);
worst->forceRadius(-1);
......
......@@ -787,7 +787,6 @@ void insertVerticesInRegion (GRegion *gr)
if(worst->isDeleted()){
myFactory.Free(worst);
allTets.erase(allTets.begin());
//Msg::Info("Worst tet is deleted");
}
else{
if(ITER++ %5000 == 0)
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment