Skip to content

Commit

Permalink
OCC: update edge/curvature meshing
Browse files Browse the repository at this point in the history
  • Loading branch information
@benoit128
benoit128 committed Dec 9, 2024
1 parent b3866b2 commit d456984
Show file tree
Hide file tree
Showing 2 changed files with 312 additions and 10 deletions.
3 changes: 2 additions & 1 deletion Cassiopee/CPlot/CPlot/Plugins/screenDump.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,13 +336,14 @@ char* Data::export2Image(E_Int exportWidth, E_Int exportHeight)
MPI_Barrier(MPI_COMM_WORLD); // seems needed

// software postprocessing on final buffer (just before screen dump)
/*
if (rank == 0)
{
char* bfl = new char [3*_view.w*_view.h];
for (E_Int i = 0; i < 3*_view.w*_view.h; i++) bfl[i] = buffer[i];
specPostProcess(bfl, _view.w, _view.h, depth, buffer);
delete [] bfl;
}
}*/

free(depth);

Expand Down
319 changes: 310 additions & 9 deletions Cassiopee/OCC/OCC/Atomic/meshEdge2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,7 @@ void geom1(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& u
//printf("h1=%f real=%f\n", h1, ue[N-1]-ue[N-2]);
}

// ============================================================================
// Geom distrib entre u0 et u1, h0 et h1/r (interieurs)
void geom2(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& ue)
{
Expand Down Expand Up @@ -186,6 +187,7 @@ void geom2(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& u
//printf("h1/r=%f real=%f\n", h1/r, ue[N-1]-ue[N-2]);
}

// ============================================================================
// Geom distrib entre u0 et u1, h0/r et h1/r (interieurs)
void geom3(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& ue)
{
Expand Down Expand Up @@ -220,6 +222,7 @@ void geom3(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& u
printf("h1/r=%f real=%f\n", h1/r, ue[N-1]-ue[N-2]);
}

// ============================================================================
// Geom distrib entre u0 et u1, h0/r et h1 (interieurs)
void geom4(E_Float u0, E_Float u1, E_Float h0, E_Float h1, E_Int& N, E_Float*& ue)
{
Expand Down Expand Up @@ -409,7 +412,8 @@ E_Int __getParamHmaxHausd(const TopoDS_Edge& E, E_Float hmax, E_Float hausd, E_I
}

// ============================================================================
// Return the nbPoints and ue for meshing E with hmin/hmax/hausd
// Return the nbPoints and ue for meshing E with hmin/hmax/hausd
// using geometric progression between points
// ============================================================================
E_Int __getParamHminHmaxHausd(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue)
{
Expand Down Expand Up @@ -551,11 +555,7 @@ E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax,

gp_Pnt point;
gp_Vec DU1, DV1, DU2, DV2, DUV;
E_Float nx, ny, nz, n2;
E_Float G11, G12, G21, G22, detG;
E_Float M11, M12, M21, M22, detM;
E_Float rho, hh, K1, K2;
E_Float det, S11, S12, S21, S22, Mi11, Mi21, Mi12, Mi22, detS;

E_Int npts = 2;
std::vector<E_Float> Us(npts);
Expand Down Expand Up @@ -586,16 +586,21 @@ E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax,
{
ues = Us[i]*(pEnd-pFirst)+pFirst;

// must be MINE (0) or OCC (1)
#define CURVATURE 1
#if CURVATURE == 0
pCurve->D0(ues, Puv);
u = Puv.X(); v = Puv.Y(); // u,v edge on surface
if (u <= BU1) u = BU1+1.e-1;
if (u >= BU2) u = BU2-1.e-1;
if (v <= BV1) v = BV1+1.e-1;
if (v >= BV2) v = BV2-1.e-1;

// must be MINE (0) or OCC (1)
#define CURVATURE 1
#if CURVATURE == 0
E_Float nx, ny, nz, n2;
E_Float G11, G12, G21, G22, detG;
E_Float M11, M12, M21, M22, detM;
E_Float det, S11, S12, S21, S22, Mi11, Mi21, Mi12, Mi22, detS;

//printf("surf: u=%g, v=%g\n", u, v);
// estimer les derivees sur la surface
surface->D2(u, v, point, DU1, DV1, DU2, DV2, DUV);
Expand Down Expand Up @@ -767,6 +772,301 @@ E_Int __getParamHminHmaxHausd2(const TopoDS_Edge& E, E_Float hmin, E_Float hmax,
return 1;

}

// ===============================================================================
// Return the nbPoints and ue for meshing E with hmin/hmax/hausd with max of faces
// new in between points evaluation
// ===============================================================================
E_Int __getParamHminHmaxHausd3(const TopoDS_Edge& E, E_Float hmin, E_Float hmax, E_Float hausd, E_Int& nbPoints, E_Float*& ue,
TopoDS_Shape* shape)
{
TopTools_IndexedDataMapOfShapeListOfShape edgeFaceMap;
TopExp::MapShapesAndAncestors(*shape, TopAbs_EDGE, TopAbs_FACE, edgeFaceMap);

BRepAdaptor_Curve C0(E);
GeomAdaptor_Curve geomAdap(C0.Curve());
Standard_Real u0 = geomAdap.FirstParameter();
Standard_Real u1 = geomAdap.LastParameter();
//printf("===============edge =========================\n");
//printf("u0=%g u1=%g\n", u0, u1);
if (BRep_Tool::Degenerated(E))
{
nbPoints = 2;
ue = new E_Float [nbPoints];
for (E_Int i = 0; i < nbPoints; i++) ue[i] = u0;
return 1;
}
E_Float L = (E_Float) GCPnts_AbscissaPoint::Length(geomAdap, u0, u1);
L = K_FUNC::E_max(L, 1.e-12);

gp_Pnt2d Puv;
E_Float u,v;
E_Float ues;
Standard_Real aFirst = C0.FirstParameter(), aEnd=C0.LastParameter();
//printf("afirst=%g aend=%g\n", aFirst, aEnd);
Standard_Real pFirst = aFirst, pEnd=aEnd;

gp_Pnt point;
gp_Vec DU1, DV1, DU2, DV2, DUV;
E_Float rho, hh, K1, K2;

E_Int npts = 3; // number of evaluation points
std::vector<E_Float> Us(npts);
for (E_Int i = 0; i < npts; i++) Us[i] = i/(npts-1.);
//Us[0] = Us[0]+0.001*(Us[1]-Us[0]);
//Us[npts-1] = Us[npts-1]-0.001*(Us[npts-1]-Us[npts-2]);

std::vector<E_Float> h(npts);
for (E_Int i = 0; i < npts; i++) h[i] = K_CONST::E_MAX_FLOAT;

// Find faces of E
const TopTools_ListOfShape& connectedFaces = edgeFaceMap.FindFromKey(E);
for (TopTools_ListIteratorOfListOfShape it(connectedFaces); it.More(); it.Next())
{
const TopoDS_Face& F = TopoDS::Face(it.Value());
Handle(Geom_Surface) surface = BRep_Tool::Surface(F);

Standard_Real BU1, BU2, BV1, BV2; surface->Bounds(BU1, BU2, BV1, BV2);
//printf("face bound = %g %g and %g %g\n", BU1, BU2, BV1, BV2);

// p curve on F
Handle(Geom2d_Curve) pCurve = BRep_Tool::CurveOnSurface(E, F, pFirst, pEnd);
//printf("pfirst=%g pend=%g\n", pFirst, pEnd);
//printf("pCurve = %p\n", (void*)pCurve);
// estimate at some points
for (E_Int i = 0; i < npts; i++)
{
ues = Us[i]*(pEnd-pFirst)+pFirst;

pCurve->D0(ues, Puv);
u = Puv.X(); v = Puv.Y(); // u,v edge on surface
if (u <= BU1) u = BU1+1.e-1;
if (u >= BU2) u = BU2-1.e-1;
if (v <= BV1) v = BV1+1.e-1;
if (v >= BV2) v = BV2-1.e-1;

// must be MINE (0) or OCC (1)
#define CURVATURE 1
#if CURVATURE == 0
E_Float nx, ny, nz, n2;
E_Float G11, G12, G21, G22, detG;
E_Float M11, M12, M21, M22, detM;
E_Float det, S11, S12, S21, S22, Mi11, Mi21, Mi12, Mi22, detS;

//printf("surf: u=%g, v=%g\n", u, v);
// estimer les derivees sur la surface
surface->D2(u, v, point, DU1, DV1, DU2, DV2, DUV);
//printf("point=%g %g %g\n", point.X(), point.Y(), point.Z());
//printf("dU1=%g %g %g\n", DU1.X(), DU1.Y(), DU1.Z());
//printf("dV1=%g %g %g\n", DV1.X(), DV1.Y(), DV1.Z());

// get normal
nx = DU1.Y()*DV1.Z()-DU1.Z()*DV1.Y();
ny = DU1.Z()*DV1.X()-DU1.X()*DV1.Z();
nz = DU1.X()*DV1.Y()-DU1.Y()*DV1.X();
n2 = nx*nx+ny*ny+nz*nz;
n2 = std::sqrt(n2);
if (n2 < 1.e-24) n2 = 1.e12;
else n2 = 1./n2;
nx = nx*n2;
ny = ny*n2;
nz = nz*n2;
//printf("normale= %g %g %g norm=%g\n", nx, ny, nz, nx*nx+ny*ny+nz*nz);

// Metric
M11 = DU1.X()*DU1.X()+DU1.Y()*DU1.Y()+DU1.Z()*DU1.Z();
M22 = DV1.X()*DV1.X()+DV1.Y()*DV1.Y()+DV1.Z()*DV1.Z();
M12 = M21 = DU1.X()*DV1.X()+DU1.Y()*DV1.Y()+DU1.Z()*DV1.Z();
//printf("matrice M:\n");
//printf(" %g %g\n", M11, M12);
//printf(" %g %g\n", M21, M22);

detM = M11*M22-M12*M21;
if (std::abs(detM) < 1.e-12) det = 1.e12; // invalid metric
else det = 1./detM;
Mi11 = det*M22;
Mi21 = -det*M21;
Mi12 = -det*M12;
Mi22 = det*M11;
detM = std::abs(detM);
//printf("detM=%g\n", detM);

// matrice de courbure
G11 = DU2.X()*nx+DU2.Y()*ny+DU2.Z()*nz;
G22 = DV2.X()*nx+DV2.Y()*ny+DV2.Z()*nz;
G12 = G21 = DUV.X()*nx+DUV.Y()*ny+DUV.Z()*nz;
detG = G11*G22-G12*G21;
detG = std::abs(detG);

S11 = G11*Mi11 + G12*Mi21;
S12 = G11*Mi12 + G12*Mi22;
S21 = G21*Mi11 + G22*Mi21;
S22 = G21*Mi12 + G22*Mi22;

detS = S11*S22-S12*S21;
//printf("dU2=%g %g %g\n", DU2.X(), DU2.Y(), DU2.Z());
//printf("dV2=%g %g %g\n", DV2.X(), DV2.Y(), DV2.Z());
//printf("dUV=%g %g %g\n", DUV.X(), DUV.Y(), DUV.Z());
//printf("detG=%g\n", detG);

// courbure de gauss - rayon K1*K2
//if (detS > 1.e-12) rho = 1./detS;
//else rho = 1.e12;
//printf("gauss=%g ou bien %g\n", detG/detM, detS);

// K1 et K2
//printf("matrice G:\n");
//printf(" %g %g\n", G11, G12);
//printf(" %g %g\n", G21, G22);

K_LINEAR::DelaunayMath sl;
sl.eigen_values(S11, S22, S12, K2, K1);
//printf("K1=%g K2=%g K1K2=%g\n", K1, K2, K1*K2);
//printf("mean=%g\n", 0.5*(S11+S22));

rho = K_FUNC::E_max(std::abs(K1), std::abs(K2));
if (rho > 1.e-12) rho = 1./rho;
else rho = 1.e12;

#else
// with open cascade
GeomLProp_SLProps props(surface, u, v, 2, Precision::Confusion());

if (props.IsCurvatureDefined())
{
Standard_Real gaussianCurvature = props.GaussianCurvature();
Standard_Real meanCurvature = props.MeanCurvature();

//std::cout << "OCC: gauss: " << gaussianCurvature << std::endl;
//std::cout << "OCC: mean: " << meanCurvature << std::endl;

// Principal curvatures
K1 = meanCurvature + sqrt(meanCurvature * meanCurvature - gaussianCurvature);
K2 = meanCurvature - sqrt(meanCurvature * meanCurvature - gaussianCurvature);
//std::cout << "OCC: k1: " << K1 << " k2: " <<K2 << std::endl;
rho = K_FUNC::E_max(std::abs(K1), std::abs(K2));
if (rho > 1.e-12) rho = 1./rho;
else rho = 1.e12;
}
else
{
rho = 1.e12;
//std::cout << "OCC: Curvature is not defined at the given parameters." << std::endl;
}
#endif

hh = std::sqrt(8.*hausd*rho);
hh = K_FUNC::E_min(hh, hmax);
hh = K_FUNC::E_max(hh, hmin);
//printf("rho=%g h=%g, h=%g\n", rho, std::sqrt(1.*hausd*rho), hh);
h[i] = K_FUNC::E_min(h[i], hh);
}
}

//printf("hi dimensioned= ");
//for (E_Int i = 0; i < npts; i++)
//{ printf("%g ", h[i]); }
//printf("\n");

for (E_Int i = 0; i < npts; i++)
{
// passage en parametres
h[i] = h[i]/L;
}

//printf("h parameter= ");
//for (E_Int i = 0; i < npts; i++)
//{ printf("%g ", h[i]); }
//printf("\n");

// compute integral of h dxsi
E_Float h1, h2, Nf, res, xi, r;
E_Int N, i0, i1;
E_Float dxi = 1./(npts-1);

E_Float integ = 0.;
for (E_Int i = 0; i < npts-1; i++)
{
h1 = h[i]; h2 = h[i+1];
integ += dxi*(h1+h2)*0.5; // mid
}
//printf("integ=%g\n", integ);

Nf = 1./integ;
N = std::rint(Nf);
N = std::max(N, 2);

//printf("N= %d, Nf = %g\n", N, Nf);

// scale h to match N
res = Nf/N;
if (std::abs(res) > 1.e-10)
{
for (E_Int i = 0; i < npts; i++) h[i] = h[i]*res;
}

// DBX
/*
integ = 0.;
for (E_Int i = 0; i < npts-1; i++)
{
h1 = h[i]; h2 = h[i+1];
integ += dxi*(h1+h2)*0.5; // mid
}
Nf = 1./integ+1.;
printf("Nf after rescale=%g\n", Nf);
*/
// END DBX

// get linear interpolation of h on new discretization
std::vector<E_Float> hn(N);
dxi = 1./(N-1);
for (E_Int j = 0; j < N; j++)
{
xi = j*dxi;
i0 = xi*(npts-1);
i1 = std::min(i0+1, npts-1);
r = xi-i0/(npts-1);
hn[j] = (1.-r)*h[i0]+r*h[i1];
}

//printf("hn remesh\n");
//for (E_Int i = 0; i < N; i++)
//{ printf("%g ", hn[i]); }
//printf("\n");

// for each point, get x
ue = new E_Float [N];

ue[0] = 0.;
for (E_Int j = 1; j < N; j++)
{
ue[j] = ue[j-1]+N*hn[j]*dxi;
}

if (std::abs(ue[N-1]-1.) > 1.e-10)
{
// rescale
res = 1./ue[N-1];
for (E_Int j = 1; j < N; j++)
{
ue[j] = ue[j]*res;
}
}

//printf("ue end\n");
//for (E_Int i = 0; i < N; i++)
//{ printf("%g ", ue[i]); }
//printf("\n");

for (E_Int j = 0; j < N; j++)
{
ue[j] = ue[j]*(aEnd-aFirst)+ aFirst;
}
nbPoints = N;
return 1;
}

// ============================================================================
// Return ue for meshing with given param in [0,1]
// ============================================================================
Expand Down Expand Up @@ -964,7 +1264,8 @@ PyObject* K_OCC::meshOneEdge(PyObject* self, PyObject* args)
else if (useFaces && hmax > 0 && hausd > 0 && externalEdge == Py_None) // mix hmax + hausd
{
// courbure sur le max des faces
__getParamHminHmaxHausd2(E, hmin, hmax, hausd, nbPoints, ue, shape);
//__getParamHminHmaxHausd2(E, hmin, hmax, hausd, nbPoints, ue, shape);
__getParamHminHmaxHausd3(E, hmin, hmax, hausd, nbPoints, ue, shape);
PyObject* o = K_ARRAY::buildArray2(4, "x,y,z,u", nbPoints, 1, 1, 1);
FldArrayF* f; K_ARRAY::getFromArray2(o, f);
__meshEdge(E, nbPoints, ue, *f, false);
Expand Down

0 comments on commit d456984

Please sign in to comment.