1 #ifndef SMOOTH_FACE_EVAL_HPP
2 #define SMOOTH_FACE_EVAL_HPP
3
4 #include "moab/Interface.hpp"
5 #include "moab/Range.hpp"
6 #include "moab/CartVect.hpp"
7 #include "MBTagConventions.hpp"
8 #include "moab/Types.hpp"
9
10 #include <cmath>
11 #include <vector>
12 #include <map>
13
14 #define determ3( p1, q1, p2, q2, p3, q3 ) \
15 ( ( q3 ) * ( ( p2 ) - ( p1 ) ) + ( q2 ) * ( ( p1 ) - ( p3 ) ) + ( q1 ) * ( ( p3 ) - ( p2 ) ) )
16 #define blend( x ) ( -2.0 * ( x ) * ( x ) * ( x ) + 3.0 * ( x ) * ( x ) )
17
18 // work only with CAMAL > = 500
19 // #if CAMAL_VERSION < 500
20 // #else
21
22 #include "moab/GeomTopoTool.hpp"
23
24 namespace moab
25 {
26 class SmoothCurve; // it is derived from SmoothBase, maybe just need
27
28 //! Implement CAMAL geometry callbacks using smooth iMesh
29 class SmoothFace // public CMLSurfEval, public SmoothBase
30 {
31 public:
32 SmoothFace( Interface* mb, EntityHandle surface_set,
33 GeomTopoTool* gTool ); // entity or entity set
34
35 virtual ~SmoothFace();
36
37 virtual double area();
38
39 virtual void bounding_box( double box_min[3], double box_max[3] );
40
41 virtual void move_to_surface( double& x, double& y, double& z );
42 /*
43 virtual void move_to_surface(double& x, double& y, double& z,
44 double& u_guess, double& v_guess);*/
45
46 virtual bool normal_at( double x, double y, double z, double& nx, double& ny, double& nz );
47
48 // initialize normals// they will be stored as tags on nodes
49 int init_gradient();
50
51 // some functions for edge evaluations
52 ErrorCode evaluate_smooth_edge( EntityHandle eh, double& t, CartVect& outv );
53
54 ErrorCode eval_bezier_patch( EntityHandle tri, CartVect& areacoord, CartVect& pt );
55
56 void project_to_facet_plane( EntityHandle tri, CartVect& pt, CartVect& point_on_plane, double& dist_to_plane );
57
58 void facet_area_coordinate( EntityHandle facet, CartVect& pt_on_plane, CartVect& areacoord );
59
60 ErrorCode project_to_facets_main( CartVect& this_point,
61 bool trim,
62 bool& outside,
63 CartVect* closest_point_ptr = NULL, // interested only in normal
64 CartVect* normal_ptr = NULL ); // interested only in closest point
65
66 ErrorCode project_to_facets( std::vector< EntityHandle >& facet_list,
67 EntityHandle& lastFacet,
68 int interpOrder,
69 double compareTol,
70 CartVect& this_point,
71 bool trim,
72 bool& outside,
73 CartVect* closest_point_ptr,
74 CartVect* normal_ptr );
75
76 ErrorCode project_to_facet( EntityHandle facet,
77 CartVect& pt,
78 CartVect& areacoord,
79 CartVect& close_point,
80 bool& outside_facet,
81 double compare_tol );
82
83 bool is_at_vertex( EntityHandle facet, // (IN) facet we are evaluating
84 CartVect& pt, // (IN) the point
85 CartVect& ac, // (IN) the ac of the point on the facet plane
86 double compare_tol, // (IN) return TRUE of closer than this
87 CartVect& eval_pt, // (OUT) location at vertex if TRUE
88 CartVect* eval_norm_ptr ); // (OUT) normal at vertex if TRUE
89
90 ErrorCode project_to_patch( EntityHandle facet, // (IN) the facet where the patch is defined
91 CartVect& ac, // (IN) area coordinate initial guess (from linear facet)
92 CartVect& pt, // (IN) point we are projecting to patch
93 CartVect& eval_pt, // (OUT) The projected point
94 CartVect* eval_norm, // (OUT) normal at evaluated point
95 bool& outside, // (OUT) the closest point on patch to pt is on an edge
96 double compare_tol, // (IN) comparison tolerance
97 int edge_id ); // (IN) only used if this is to be projected to one
98 // of the edges. Otherwise, should be -1
99
100 ErrorCode eval_bezier_patch_normal( EntityHandle facet, CartVect& areacoord, CartVect& normal );
101
102 // this will be called now from driver...
103 ErrorCode compute_tangents_for_each_edge(); // they will be used for control points
104
105 ErrorCode get_normals_for_vertices( const EntityHandle* conn2,
106 CartVect N[2] ); // here we need the gradient tag
107
108 // make this one public, will be called during initialization !!!
109 ErrorCode init_edge_control_points( CartVect& P0,
110 CartVect& P3,
111 CartVect& N0,
112 CartVect& N3,
113 CartVect& T0,
114 CartVect& T3,
115 CartVect* Pi );
116
117 // moved from private, because will be called from PaveDriver
118 ErrorCode compute_control_points_on_edges( double min_dot, Tag edgeCtrlTag, Tag markTag );
119
120 ErrorCode compute_internal_control_points_on_facets( double min_dot, Tag facetCtrlTag, Tag facetEdgeCtrlTag );
121
122 // move from private too
123 void DumpModelControlPoints();
124
125 int eval_counter()
126 {
127 return _evaluationsCounter;
128 }
129
130 // new method for ray intersection correction
131 ErrorCode ray_intersection_correct( EntityHandle facet, // (IN) the facet where the patch is defined
132 CartVect& pt, // (IN) shoot from
133 CartVect& ray, // (IN) ray direction
134 CartVect& eval_pt, // (INOUT) The intersection point
135 double& distance, // (IN OUT) the new distance
136 bool& outside ); // (OUT) the closest point on patch to pt is on an edge
137
138 void append_smooth_tags( std::vector< Tag >& smoothTags );
139 // of the edges. Otherwise, should be -1)
140 private:
141 //===========================================================================
142 // Function Name: move_ac_inside
143 //
144 // Member Type: PRIVATE
145 // Description: find the closest area coordinate to the boundary of the
146 // patch if any of its components are < 0
147 // Return if the ac was modified.
148 //===========================================================================
149 bool move_ac_inside( CartVect& ac, double tol );
150
151 //===========================================================================
152 // Function Name: ac_at_edge
153 //
154 // Member Type: PRIVATE
155 // Description: determine the area coordinate of the facet at the edge
156 //===========================================================================
157 void ac_at_edge( CartVect& fac, // facet area coordinate
158 CartVect& eac, // edge area coordinate
159 int edge_id ); // id of edge
160
161 // some local functions that do not need to be public
162 ErrorCode init_bezier_edge( EntityHandle edge, double min_dot );
163 //
164
165 ErrorCode find_edges_orientations( EntityHandle edges[3],
166 const EntityHandle* conn3,
167 int orient[3] ); // maybe we will set it?
168
169 ErrorCode init_facet_control_points( CartVect N[6], // vertex normals (per edge)
170 CartVect P[3][5], // edge control points
171 CartVect G[6] ); // return internal control points
172
173 // those are the bounding box limits;
174 // they are adjusted for the control points too in each triangle
175 void adjust_bounding_box( CartVect& vect );
176 CartVect _minim;
177 CartVect _maxim;
178
179 Range _triangles;
180 Range _edges;
181 Range _nodes;
182
183 // std::vector<double> _fractions;// they are increasing from 0. to 1., do we need these?
184 // std::vector<double> _loopLengths;
185
186 // number of loops is decided by the size of _loopEnds.size()
187 // this ref face will be gone, we will replace it with a new call
188 // RefFace * _smooth_face;
189 // int myDimension;
190 // double meshSize;
191
192 // this tag is on edges
193 // rval = _mb->tag_create("MARKER", 1, MB_TAG_BIT, _markTag, &value);
194 Tag _markTag; // this is a tag used to mark edges when we look for loops
195
196 // this tag is on nodes
197 // ErrorCode rval = _mb->tag_create("GRADIENT", 3 * sizeof(double),
198 // MB_TAG_DENSE, _gradientTag, &defNormal);
199 Tag _gradientTag; // this will be used for normal at nodes
200
201 // this tag is on edges
202 // ErrorCode rval = _mb->tag_create("TANGENTS", 6 * sizeof(double),
203 // MB_TAG_DENSE, _tangentsTag, &defTangents);
204 Tag _tangentsTag; // each edge will have exactly 2 tangents, because there is only
205 // one feature edge, and it is periodic
206 // the feature edge is exactly on the boundary
207
208 // this tag is on edges
209 // ErrorCode rval = _mb->tag_create("CONTROLEDGE", 9 * sizeof(double),
210 // MB_TAG_DENSE, _edgeCtrlTag, &defControls);
211 Tag _edgeCtrlTag;
212
213 // this tag is on facets (triangles), 6 control points on each facet
214 // there are also some 15 points used in evaluation; how to store them?
215 // ErrorCode rval = _mb->tag_create("CONTROLFACE", 18 * sizeof(double),
216 // MB_TAG_DENSE, _facetCtrlTag, &defControls);
217 Tag _facetCtrlTag;
218
219 // these are the 12 points stored for each edge, again
220 // it is cheaper this way compared to retrieve the edges every time, determine their
221 // orientation, order in triangle, and retrieve the control points from the edge the control
222 // points are stored as 12 points, in order : edge 0, 1, and 2, in that order
223 // ErrorCode rval = _mb->tag_create("CONTROLEDGEFACE", 36 * sizeof(double),
224 // MB_TAG_DENSE, _facetEdgeCtrlTag, &defControls);
225 Tag _facetEdgeCtrlTag; //
226 // plane of the facet stores as a normal a, b, c and d, distance, for
227 // ax+by+cz+d=0
228 // ErrorCode rval = _mb->tag_create("PLANE", 4 * sizeof(double),
229 // MB_TAG_DENSE, _planeTag, &defPlane);
230 Tag _planeTag;
231
232 Interface* _mb;
233 EntityHandle _set;
234 GeomTopoTool* _my_geomTopoTool;
235 EntityHandle _obb_root;
236 // counter for calls
237 long _evaluationsCounter;
238 };
239 // #endif
240 } // namespace moab
241 #endif /* SMOOTH_FACE_EVAL_HPP*/