diff --git a/Mesh/Field.cpp b/Mesh/Field.cpp index 7f937607c2a555e5e42f46c04677fac33b54c504..6888cc6a3f875d869ca8fc5c725946615642afa6 100644 --- a/Mesh/Field.cpp +++ b/Mesh/Field.cpp @@ -1330,6 +1330,140 @@ struct AttractorInfo{ double u,v; }; +class AttractorAnisoCurveField : public Field { + ANNkd_tree *kdtree; + ANNpointArray zeronodes; + ANNidxArray index; + ANNdistArray dist; + std::list<int> edges_id; + double dMin, dMax, lMinTangent, lMaxTangent, lMinNormal, lMaxNormal; + int n_nodes_by_edge; + std::vector<SVector3> tg; + public: + AttractorAnisoCurveField() : kdtree(0), zeronodes(0) + { + index = new ANNidx[1]; + dist = new ANNdist[1]; + n_nodes_by_edge = 20; + update_needed = true; + dMin = 0.1; + dMax = 0.5; + lMinNormal = 0.05; + lMinTangent = 0.5; + lMaxNormal = 0.5; + lMaxTangent = 0.5; + options["EdgesList"] = new FieldOptionList + (edges_id, "Indices of curves in the geometric model", &update_needed); + options["NNodesByEdge"] = new FieldOptionInt + (n_nodes_by_edge, "Number of nodes used to discretized each curve", + &update_needed); + options["dMin"] = new FieldOptionDouble + (dMin, "Minimum distance, bellow this distance from the curves, prescribe the minimum mesh sizes."); + options["dMax"] = new FieldOptionDouble + (dMax, "Maxmium distance, above this distance from the curves, prescribe the maximum mesh sizes."); + options["lMinTangent"] = new FieldOptionDouble + (lMinTangent, "Minimum mesh size in the direction tangeant to the closest curve."); + options["lMaxTangent"] = new FieldOptionDouble + (lMaxTangent, "Maximum mesh size in the direction tangeant to the closest curve."); + options["lMinNormal"] = new FieldOptionDouble + (lMinNormal, "Minimum mesh size in the direction normal to the closest curve."); + options["lMaxNormal"] = new FieldOptionDouble + (lMaxNormal, "Maximum mesh size in the direction normal to the closest curve."); + } + virtual bool isotropic () const {return false;} + ~AttractorAnisoCurveField() + { + if(kdtree) delete kdtree; + if(zeronodes) annDeallocPts(zeronodes); + delete[]index; + delete[]dist; + } + const char *getName() + { + return "AttractorAnisoCurve"; + } + std::string getDescription() + { + return "Compute the distance from the nearest curve in a list. Then the mesh " + "size can be specified independently in the direction normal to the curve " + "and in the direction parallel to the curve (Each curve " + "is replaced by NNodesByEdge equidistant nodes and the distance from those " + "nodes is computed.)"; + } + void update() { + if(zeronodes) { + annDeallocPts(zeronodes); + delete kdtree; + } + int totpoints = n_nodes_by_edge * edges_id.size(); + if(totpoints){ + zeronodes = annAllocPts(totpoints, 3); + } + tg.resize(totpoints); + int k = 0; + for(std::list<int>::iterator it = edges_id.begin(); + it != edges_id.end(); ++it) { + Curve *c = FindCurve(*it); + if(c) { + for(int i = 0; i < n_nodes_by_edge; i++) { + double u = (double)i / (n_nodes_by_edge - 1); + Vertex V = InterpolateCurve(c, u, 0); + zeronodes[k][0] = V.Pos.X; + zeronodes[k][1] = V.Pos.Y; + zeronodes[k][2] = V.Pos.Z; + Vertex V2 = InterpolateCurve(c, u, 1); + tg[k] = SVector3(V2.Pos.X, V2.Pos.Y, V2.Pos.Z); + tg[k].normalize(); + k++; + } + } + else { + GEdge *e = GModel::current()->getEdgeByTag(*it); + if(e) { + for(int i = 0; i < n_nodes_by_edge; i++) { + double u = (double)i / (n_nodes_by_edge - 1); + Range<double> b = e->parBounds(0); + double t = b.low() + u * (b.high() - b.low()); + GPoint gp = e->point(t); + SVector3 d = e->firstDer(t); + zeronodes[k][0] = gp.x(); + zeronodes[k][1] = gp.y(); + zeronodes[k][2] = gp.z(); + tg[k] = d; + tg[k].normalize(); + k++; + } + } + } + } + kdtree = new ANNkd_tree(zeronodes, totpoints, 3); + update_needed = false; + } + void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0) + { + if(update_needed) + update(); + double xyz[3] = { x, y, z }; + kdtree->annkSearch(xyz, 1, index, dist); + double d = sqrt(dist[0]); + double lTg = d < dMin ? lMinTangent : d > dMax ? lMaxTangent : lMinTangent + (lMaxTangent - lMinTangent) * (d - dMin) / (dMax - dMin); + double lN = d < dMin ? lMinNormal : d > dMax ? lMaxNormal : lMinNormal + (lMaxNormal - lMinNormal) * (d - dMin) / (dMax - dMin); + SVector3 t = tg[index[0]]; + SVector3 n0 = crossprod(t, fabs(t(0)) > fabs(t(1)) ? SVector3(0,1,0) : SVector3(1,0,0)); + SVector3 n1 = crossprod(t, n0); + metr = SMetric3(1/lTg/lTg, 1/lN/lN, 1/lN/lN, t, n0, n1); + } + virtual double operator() (double X, double Y, double Z, GEntity *ge=0) + { + if(update_needed) + update(); + double xyz[3] = { X, Y, Z }; + kdtree->annkSearch(xyz, 1, index, dist); + double d = sqrt(dist[0]); + return std::max(d, 0.05); + } +}; + class AttractorField : public Field { ANNkd_tree *kdtree; @@ -1689,6 +1823,7 @@ FieldManager::FieldManager() map_type_name["MathEvalAniso"] = new FieldFactoryT<MathEvalFieldAniso>(); #if defined(HAVE_ANN) map_type_name["Attractor"] = new FieldFactoryT<AttractorField>(); + map_type_name["AttractorAnisoCurve"] = new FieldFactoryT<AttractorAnisoCurveField>(); #endif map_type_name["MaxEigenHessian"] = new FieldFactoryT<MaxEigenHessianField>(); background_field = -1;