moab
|
00001 #ifndef SMOOTH_FACE_EVAL_HPP 00002 #define SMOOTH_FACE_EVAL_HPP 00003 00004 // do we really need iMesh here; maybe go directly to MOAB 00005 //#include "iMesh.h" 00006 #include "moab/Interface.hpp" 00007 #include "moab/Range.hpp" 00008 #include "moab/CartVect.hpp" 00009 #include "MBTagConventions.hpp" 00010 #include "moab/Types.hpp" 00011 00012 #define determ3(p1,q1,p2,q2,p3,q3) ((q3)*((p2)-(p1)) + (q2)*((p1)-(p3)) + (q1)*((p3)-(p2))) 00013 #define sqr(a) ((a)*(a)) 00014 #define cube(a) (sqr(a) * (a)) 00015 #define quart(a) (sqr(a) * sqr(a)) 00016 #define blend(x) (-2.0*(x)*(x)*(x) + 3.0*(x)*(x)) 00017 00018 #include <vector> 00019 #include <map> 00020 //#include "MBEntityHandle.hpp" 00021 00022 // work only with CAMAL > = 500 00023 // #if CAMAL_VERSION < 500 00024 00025 // #else 00026 00027 #include "moab/GeomTopoTool.hpp" 00028 00029 namespace moab { 00030 class SmoothCurve;// it is derived from SmoothBase, maybe just need 00031 00033 class SmoothFace // public CMLSurfEval, public SmoothBase 00034 { 00035 public: 00036 SmoothFace(Interface * mb, EntityHandle surface_set, GeomTopoTool * gTool); // entity or entity set 00037 00038 virtual ~SmoothFace(); 00039 00040 virtual double area(); 00041 00042 virtual void bounding_box(double box_min[3], double box_max[3]); 00043 00044 virtual void move_to_surface(double& x, double& y, double& z); 00045 /* 00046 virtual void move_to_surface(double& x, double& y, double& z, 00047 double& u_guess, double& v_guess);*/ 00048 00049 virtual bool normal_at(double x, double y, double z, double& nx, double& ny, 00050 double& nz); 00051 00052 // initialize normals// they will be stored as tags on nodes 00053 int init_gradient(); 00054 00055 // some functions for edge evaluations 00056 ErrorCode evaluate_smooth_edge(EntityHandle eh, double &t, CartVect & outv); 00057 00058 ErrorCode eval_bezier_patch(EntityHandle tri, CartVect &areacoord, 00059 CartVect &pt); 00060 00061 void project_to_facet_plane(EntityHandle tri, CartVect &pt, 00062 CartVect &point_on_plane, double &dist_to_plane); 00063 00064 void facet_area_coordinate(EntityHandle facet, CartVect & pt_on_plane, 00065 CartVect & areacoord); 00066 00067 ErrorCode project_to_facets_main(CartVect &this_point, bool trim, 00068 bool & outside, CartVect * closest_point_ptr = NULL, // interested only in normal 00069 CartVect * normal_ptr = NULL); // interested only in closest point 00070 00071 ErrorCode project_to_facets(std::vector<EntityHandle> & facet_list, 00072 EntityHandle & lastFacet, int interpOrder, double compareTol, 00073 CartVect &this_point, bool trim, bool & outside, 00074 CartVect *closest_point_ptr, CartVect * normal_ptr); 00075 00076 ErrorCode project_to_facet(EntityHandle facet, CartVect &pt, 00077 CartVect &areacoord, CartVect &close_point, bool &outside_facet, 00078 double compare_tol); 00079 00080 bool is_at_vertex(EntityHandle facet, // (IN) facet we are evaluating 00081 CartVect &pt, // (IN) the point 00082 CartVect &ac, // (IN) the ac of the point on the facet plane 00083 double compare_tol, // (IN) return TRUE of closer than this 00084 CartVect &eval_pt, // (OUT) location at vertex if TRUE 00085 CartVect *eval_norm_ptr); // (OUT) normal at vertex if TRUE 00086 00087 ErrorCode project_to_patch(EntityHandle facet, // (IN) the facet where the patch is defined 00088 CartVect &ac, // (IN) area coordinate initial guess (from linear facet) 00089 CartVect &pt, // (IN) point we are projecting to patch 00090 CartVect &eval_pt, // (OUT) The projected point 00091 CartVect *eval_norm, // (OUT) normal at evaluated point 00092 bool &outside, // (OUT) the closest point on patch to pt is on an edge 00093 double compare_tol, // (IN) comparison tolerance 00094 int edge_id); // (IN) only used if this is to be projected to one 00095 // of the edges. Otherwise, should be -1 00096 00097 ErrorCode eval_bezier_patch_normal(EntityHandle facet, CartVect &areacoord, 00098 CartVect &normal); 00099 00100 // this will be called now from driver... 00101 ErrorCode compute_tangents_for_each_edge();// they will be used for control points 00102 00103 ErrorCode get_normals_for_vertices(const EntityHandle * conn2, CartVect N[2]);// here we need the gradient tag 00104 00105 // make this one public, will be called during initialization !!! 00106 ErrorCode init_edge_control_points(CartVect &P0, CartVect &P3, CartVect &N0, 00107 CartVect &N3, CartVect &T0, CartVect &T3, CartVect * Pi); 00108 00109 // moved from private, because will be called from PaveDriver 00110 ErrorCode compute_control_points_on_edges(double min_dot, Tag edgeCtrlTag, 00111 Tag markTag); 00112 00113 ErrorCode compute_internal_control_points_on_facets(double min_dot, 00114 Tag facetCtrlTag, Tag facetEdgeCtrlTag); 00115 00116 // move from private too 00117 void DumpModelControlPoints(); 00118 00119 int eval_counter() 00120 { 00121 return _evaluationsCounter; 00122 } 00123 00124 // new method for ray intersection correction 00125 ErrorCode ray_intersection_correct(EntityHandle facet, // (IN) the facet where the patch is defined 00126 CartVect &pt, // (IN) shoot from 00127 CartVect &ray, // (IN) ray direction 00128 CartVect &eval_pt, // (INOUT) The intersection point 00129 double & distance, // (IN OUT) the new distance 00130 bool &outside); // (OUT) the closest point on patch to pt is on an edge 00131 00132 void append_smooth_tags(std::vector<Tag> & smoothTags); 00133 // of the edges. Otherwise, should be -1) 00134 private: 00135 00136 //=========================================================================== 00137 //Function Name: move_ac_inside 00138 // 00139 //Member Type: PRIVATE 00140 //Description: find the closest area coordinate to the boundary of the 00141 // patch if any of its components are < 0 00142 // Return if the ac was modified. 00143 //=========================================================================== 00144 bool move_ac_inside(CartVect &ac, double tol); 00145 00146 //=========================================================================== 00147 //Function Name: ac_at_edge 00148 // 00149 //Member Type: PRIVATE 00150 //Description: determine the area coordinate of the facet at the edge 00151 //=========================================================================== 00152 void ac_at_edge(CartVect &fac, // facet area coordinate 00153 CartVect &eac, // edge area coordinate 00154 int edge_id); // id of edge 00155 00156 // some local functions that do not need to be public 00157 ErrorCode init_bezier_edge(EntityHandle edge, double min_dot); 00158 // 00159 00160 ErrorCode find_edges_orientations(EntityHandle edges[3], 00161 const EntityHandle * conn3, int orient[3]); // maybe we will set it? 00162 00163 ErrorCode init_facet_control_points(CartVect N[6], // vertex normals (per edge) 00164 CartVect P[3][5], // edge control points 00165 CartVect G[6]); // return internal control points 00166 00167 // those are the bounding box limits; 00168 // they are adjusted for the control points too in each triangle 00169 void adjust_bounding_box(CartVect & vect); 00170 CartVect _minim; 00171 CartVect _maxim; 00172 00173 Range _triangles; 00174 Range _edges; 00175 Range _nodes; 00176 00177 //std::vector<double> _fractions;// they are increasing from 0. to 1., do we need these? 00178 //std::vector<double> _loopLengths; 00179 00180 // number of loops is decided by the size of _loopEnds.size() 00181 // this ref face will be gone, we will replace it with a new call 00182 //RefFace * _smooth_face; 00183 //int myDimension; 00184 //double meshSize; 00185 00186 // this tag is on edges 00187 // rval = _mb->tag_create("MARKER", 1, MB_TAG_BIT, _markTag, &value); 00188 Tag _markTag; // this is a tag used to mark edges when we look for loops 00189 00190 // this tag is on nodes 00191 //ErrorCode rval = _mb->tag_create("GRADIENT", 3 * sizeof(double), 00192 // MB_TAG_DENSE, _gradientTag, &defNormal); 00193 Tag _gradientTag; // this will be used for normal at nodes 00194 00195 // this tag is on edges 00196 //ErrorCode rval = _mb->tag_create("TANGENTS", 6 * sizeof(double), 00197 // MB_TAG_DENSE, _tangentsTag, &defTangents); 00198 Tag _tangentsTag; // each edge will have exactly 2 tangents, because there is only 00199 // one feature edge, and it is periodic 00200 // the feature edge is exactly on the boundary 00201 00202 // this tag is on edges 00203 //ErrorCode rval = _mb->tag_create("CONTROLEDGE", 9 * sizeof(double), 00204 // MB_TAG_DENSE, _edgeCtrlTag, &defControls); 00205 Tag _edgeCtrlTag; 00206 00207 // this tag is on facets (triangles), 6 control points on each facet 00208 // there are also some 15 points used in evaluation; how to store them? 00209 //ErrorCode rval = _mb->tag_create("CONTROLFACE", 18 * sizeof(double), 00210 // MB_TAG_DENSE, _facetCtrlTag, &defControls); 00211 Tag _facetCtrlTag; 00212 00213 // these are the 12 points stored for each edge, again 00214 // it is cheaper this way compared to retrieve the edges every time, determine their orientation, order 00215 // in triangle, and retrieve the control points from the edge 00216 // the control points are stored as 12 points, in order : edge 0, 1, and 2, in that order 00217 //ErrorCode rval = _mb->tag_create("CONTROLEDGEFACE", 36 * sizeof(double), 00218 // MB_TAG_DENSE, _facetEdgeCtrlTag, &defControls); 00219 Tag _facetEdgeCtrlTag; // 00220 // plane of the facet stores as a normal a, b, c and d, distance, for 00221 // ax+by+cz+d=0 00222 //ErrorCode rval = _mb->tag_create("PLANE", 4 * sizeof(double), 00223 // MB_TAG_DENSE, _planeTag, &defPlane); 00224 Tag _planeTag; 00225 00226 Interface * _mb; 00227 EntityHandle _set; 00228 GeomTopoTool * _my_geomTopoTool; 00229 EntityHandle _obb_root; 00230 // counter for calls 00231 long _evaluationsCounter; 00232 }; 00233 // #endif 00234 } // namespace moab 00235 #endif /* SMOOTH_FACE_EVAL_HPP*/