Branch data Line data Source code
1 : : //
2 : : // File: CubitQuadFacet.cpp
3 : : //
4 : : // Owner: sjowen
5 : : //
6 : : #include "CubitDefines.h"
7 : : #include "CubitQuadFacet.hpp"
8 : : #include "CubitFacetEdge.hpp"
9 : : #include "CubitPoint.hpp"
10 : : #include "CubitFacet.hpp"
11 : : #include "CubitVector.hpp"
12 : : #include "FacetEvalTool.hpp"
13 : :
14 : : #define DETERM(p1,q1,p2,q2,p3,q3) ((q3)*((p2)-(p1)) + (q2)*((p1)-(p3)) + (q1)*((p3)-(p2)))
15 : :
16 : : //===========================================================================
17 : : // Function: CubitQuadFacet
18 : : // Purpose: constructor
19 : : // Date: 4/11/01
20 : : // Author: sjowen
21 : : //===========================================================================
22 : 0 : CubitQuadFacet::CubitQuadFacet()
23 : : {
24 : 0 : }
25 : :
26 : : //===========================================================================
27 : : // Function: CubitQuadFacet
28 : : // Purpose: destructor
29 : : // Date: 4/11/01
30 : : // Author: sjowen
31 : : //===========================================================================
32 : 0 : CubitQuadFacet::~CubitQuadFacet()
33 : : {
34 [ # # ]: 0 : }
35 : :
36 : : //===========================================================================
37 : : // Function: evaluate
38 : : // Purpose: evaluate the quad at a specified u,v coordinate
39 : : // Notes: Uses the interpolation order previously set up for the
40 : : // facet-based surface. The u,v system is based on the point
41 : : // order and is defined as follows:
42 : : //
43 : : // [3] = (-1,1) *============* [2] = (1,1)
44 : : // | |
45 : : // | |
46 : : // | |
47 : : // | |
48 : : // | |
49 : : // [0] = (-1,-1) *============* [1] = (1,-1)
50 : : //
51 : : // Date: 4/11/01
52 : : // Author: sjowen
53 : : //===========================================================================
54 : 0 : CubitStatus CubitQuadFacet::evaluate( double u, double v,
55 : : CubitVector *eval_point,
56 : : CubitVector *eval_normal,
57 : : CubitVector *eval_du,
58 : : CubitVector *eval_dv)
59 : : {
60 [ # # ]: 0 : CubitVector normal;
61 : 0 : CubitStatus stat = CUBIT_SUCCESS;
62 : : int tri_index;
63 : : double A, B, C;
64 [ # # ]: 0 : map_to_tri_system( u, v, tri_index, A, B, C );
65 [ # # ]: 0 : CubitVector ac(A,B,C);
66 [ # # ]: 0 : CubitFacet *tri_facet = get_tri_facet( tri_index );
67 [ # # ]: 0 : CubitVector point;
68 [ # # ]: 0 : stat = tri_facet->evaluate( ac, &point, eval_normal );
69 [ # # ]: 0 : if (!stat)
70 : 0 : return stat;
71 [ # # ]: 0 : if (eval_point != NULL)
72 [ # # ]: 0 : *eval_point = point;
73 [ # # ][ # # ]: 0 : if (eval_du != NULL || eval_dv != NULL)
74 : : {
75 [ # # ][ # # ]: 0 : CubitVector du, dv;
76 [ # # ]: 0 : stat = eval_derivatives( u, v, point, du, dv );
77 [ # # ][ # # ]: 0 : normal = du * dv;
78 [ # # ]: 0 : normal.normalize();
79 [ # # ]: 0 : if (eval_du != NULL)
80 [ # # ]: 0 : *eval_du = du;
81 [ # # ]: 0 : if (eval_dv != NULL)
82 [ # # ]: 0 : *eval_dv = dv;
83 [ # # ]: 0 : if (eval_normal != NULL)
84 [ # # ]: 0 : *eval_normal = normal;
85 : : }
86 : 0 : return stat;
87 : : }
88 : :
89 : : //===========================================================================
90 : : // Function: eval_derivatives
91 : : // Purpose: evaluate derivatives on quad with respect to the u,v system
92 : : // on the quad
93 : : // Date: 5/01
94 : : // Author: sjowen
95 : : //===========================================================================
96 : 0 : CubitStatus CubitQuadFacet::eval_derivatives( double u, double v,
97 : : CubitVector &point,
98 : : CubitVector &du,
99 : : CubitVector &dv )
100 : : {
101 : : // do a forward, backward or central difference depending on where the
102 : : // u, v is located on the quad
103 : :
104 : 0 : CubitStatus stat = CUBIT_SUCCESS;
105 : 0 : double incr = 0.001;
106 : 0 : double min = -1.0;
107 : 0 : double max = 1.0;
108 : : double uf, vf, ub, vb;
109 [ # # ][ # # ]: 0 : CubitVector pf, pb;
110 [ # # ]: 0 : if (u < min + incr)
111 : : {
112 : 0 : uf = u + incr;
113 [ # # ]: 0 : stat = evaluate(uf, v, &pf);
114 [ # # ][ # # ]: 0 : du = ( pf - point ) / incr;
[ # # ]
115 : : }
116 [ # # ]: 0 : else if (u > max - incr)
117 : : {
118 : 0 : ub = u - incr;
119 [ # # ]: 0 : stat = evaluate( ub, v, &pb);
120 [ # # ][ # # ]: 0 : du = ( point - pb ) / incr;
[ # # ]
121 : : }
122 : : else
123 : : {
124 : 0 : uf = u + incr;
125 : 0 : ub = u - incr;
126 [ # # ]: 0 : stat = evaluate( uf, v, &pf );
127 [ # # ]: 0 : stat = evaluate( ub, v, &pb );
128 [ # # ][ # # ]: 0 : du = ( pf - pb ) / (2.0 * incr);
[ # # ]
129 : : }
130 [ # # ]: 0 : if (v < min + incr)
131 : : {
132 : 0 : vf = v + incr;
133 [ # # ]: 0 : stat = evaluate(u, vf, &pf);
134 [ # # ][ # # ]: 0 : dv = ( pf - point ) / incr;
[ # # ]
135 : : }
136 [ # # ]: 0 : else if (v > max - incr)
137 : : {
138 : 0 : vb = v - incr;
139 [ # # ]: 0 : stat = evaluate( u, vb, &pb);
140 [ # # ][ # # ]: 0 : dv = ( point - pb ) / incr;
[ # # ]
141 : : }
142 : : else
143 : : {
144 : 0 : vf = v + incr;
145 : 0 : vb = v - incr;
146 [ # # ]: 0 : stat = evaluate( u, vf, &pf );
147 [ # # ]: 0 : stat = evaluate( u, vb, &pb );
148 [ # # ][ # # ]: 0 : dv = ( pf - pb ) / (2.0 * incr);
[ # # ]
149 : : }
150 : 0 : return stat;
151 : : }
152 : :
153 : : //===========================================================================
154 : : // Function: debug_draw
155 : : // Purpose: draw the facets
156 : : // Date: 4/11/01
157 : : // Author: sjowen
158 : : //===========================================================================
159 : 0 : void CubitQuadFacet::debug_draw( int color, int flush_it, int draw_uv )
160 : : {
161 : 0 : CubitFacet *tri_facet = get_tri_facet( 0 );
162 : 0 : tri_facet->debug_draw(color, flush_it, draw_uv);
163 : 0 : tri_facet = get_tri_facet( 1 );
164 : 0 : tri_facet->debug_draw(color, flush_it, draw_uv);
165 : 0 : }
166 : :
167 : : //===========================================================================
168 : : // Function: map_to_tri_system
169 : : // Purpose: given a u,v coordinate, return the triangle and area coordinate
170 : : // Date: 4/11/01
171 : : // Author: sjowen
172 : : //===========================================================================
173 : 0 : void CubitQuadFacet::map_to_tri_system(double u, double v,
174 : : int &tri_index,
175 : : double &A, double &B, double &C)
176 : : {
177 : : double trivert[3][2];
178 : : double areacoords[2][3];
179 : : int ii, jj, j0, j1;
180 : : double small_ac[2];
181 : 0 : small_ac[0] = small_ac[1] = 1.0;
182 : :
183 : : // loop through both tris
184 : :
185 [ # # ]: 0 : for (ii=0; ii<2; ii++)
186 : : {
187 : :
188 : : // set up the (u,v) vertex coordinates for the tri
189 : :
190 [ # # ]: 0 : for (jj=0; jj<3; jj++)
191 : : {
192 [ # # ][ # # : 0 : switch(tri_to_quad_index(ii,jj))
# # # ]
193 : : {
194 : : case 0:
195 : 0 : trivert[jj][0] = -1.0;
196 : 0 : trivert[jj][1] = -1.0;
197 : 0 : break;
198 : : case 1:
199 : 0 : trivert[jj][0] = 1.0;
200 : 0 : trivert[jj][1] = -1.0;
201 : 0 : break;
202 : : case 2:
203 : 0 : trivert[jj][0] = 1.0;
204 : 0 : trivert[jj][1] = 1.0;
205 : 0 : break;
206 : : case 3:
207 : 0 : trivert[jj][0] = -1.0;
208 : 0 : trivert[jj][1] = 1.0;
209 : 0 : break;
210 : : }
211 : : }
212 : :
213 : : // compute the area coordinates
214 : :
215 [ # # ]: 0 : for (jj=0; jj<3; jj++)
216 : : {
217 : 0 : j0 = (jj+1)%3;
218 : 0 : j1 = (jj+2)%3;
219 : 0 : areacoords[ii][jj] = DETERM( trivert[j0][0], trivert[j0][1],
220 : 0 : trivert[j1][0], trivert[j1][1], u, v );
221 [ # # ]: 0 : if (small_ac[ii] > areacoords[ii][jj])
222 : : {
223 : 0 : small_ac[ii] = areacoords[ii][jj];
224 : : }
225 : : }
226 : : }
227 : 0 : tri_index = (small_ac[0] > small_ac[1]) ? 0 : 1;
228 : 0 : A = areacoords[tri_index][0];
229 : 0 : B = areacoords[tri_index][1];
230 : 0 : C = areacoords[tri_index][2];
231 : 0 : double sum = A + B + C;
232 : 0 : A /= sum;
233 : 0 : B /= sum;
234 : 0 : C /= sum;
235 : 0 : }
236 : :
237 : : //===========================================================================
238 : : // Function: get_control_points
239 : : // Purpose: retreive the control points from the quad
240 : : // Date: 4/11/01
241 : : // Author: sjowen
242 : : //===========================================================================
243 : 0 : void CubitQuadFacet::get_control_points( double *ctrl_pts )
244 : : {
245 : : // assumes 6 control points per triangle and 3 along the diagonal edge
246 : :
247 : : int ii, jj;
248 : 0 : int ipt = 0;
249 : 0 : int index = 0;
250 : 0 : int numtri = 2;
251 : 0 : int num_ctrl_pts_per_tri = 6;
252 [ # # ]: 0 : for (ii=0; ii<numtri; ii++)
253 : : {
254 : 0 : CubitFacet *tri_facet = get_tri_facet( ii );
255 : 0 : CubitVector *tri_ctrl_pts = tri_facet->control_points();
256 [ # # ]: 0 : assert(tri_ctrl_pts != NULL);
257 [ # # ]: 0 : for (jj=0; jj<num_ctrl_pts_per_tri; jj++)
258 : : {
259 : 0 : index = ipt * 3;
260 : 0 : ctrl_pts[index] = tri_ctrl_pts[jj].x();
261 : 0 : ctrl_pts[index+1] = tri_ctrl_pts[jj].y();
262 : 0 : ctrl_pts[index+2] = tri_ctrl_pts[jj].z();
263 : 0 : ipt++;
264 : : }
265 : : }
266 : :
267 : 0 : int num_ctrl_pts_per_edge = 3;
268 : 0 : CubitVector *edge_ctrl_pts = NULL;
269 : 0 : CubitFacetEdge *edge_ptr = this->point(0)->shared_edge( this->point(2) );
270 [ # # ]: 0 : assert(edge_ptr != NULL);
271 : 0 : edge_ctrl_pts = edge_ptr->control_points();
272 [ # # ]: 0 : for (ii=0; ii<num_ctrl_pts_per_edge; ii++)
273 : : {
274 : 0 : index = ipt * 3;
275 : 0 : ctrl_pts[index] = edge_ctrl_pts[ii].x();
276 : 0 : ctrl_pts[index+1] = edge_ctrl_pts[ii].y();
277 : 0 : ctrl_pts[index+2] = edge_ctrl_pts[ii].z();
278 : 0 : ipt++;
279 : : }
280 : 0 : }
281 : :
282 : : //===========================================================================
283 : : // Function: set_control_points
284 : : // Purpose: set the control points into the quad
285 : : // Date: 4/11/01
286 : : // Author: sjowen
287 : : //===========================================================================
288 : 0 : void CubitQuadFacet::set_control_points( double *ctrl_pts )
289 : : {
290 : : // assumes 6 control points per triangle and 3 along the diagonal edge
291 : :
292 : : int ii;
293 : 0 : int index = 0;
294 : 0 : int numtri = 2;
295 : 0 : int num_ctrl_pts_per_tri = 6;
296 [ # # ]: 0 : for (ii=0; ii<numtri; ii++)
297 : : {
298 : 0 : CubitFacet *tri_facet = get_tri_facet( ii );
299 : 0 : index = 3 * ii * num_ctrl_pts_per_tri;
300 : 0 : tri_facet->set_control_points( &(ctrl_pts[index]) );
301 : : }
302 : :
303 : 0 : CubitFacetEdge *edge_ptr = this->point(0)->shared_edge( this->point(2) );
304 [ # # ]: 0 : assert(edge_ptr != NULL);
305 : 0 : index = 3 * numtri * num_ctrl_pts_per_tri;
306 : 0 : edge_ptr->set_control_points( &(ctrl_pts[index]) );
307 : 0 : }
308 : :
309 : : //===========================================================================
310 : : // Function: points
311 : : // Purpose: get the points on the quad. Virtual to FacetEntity
312 : : // Date: 8/1/03
313 : : // Author: sjowen
314 : : //===========================================================================
315 : 0 : void CubitQuadFacet::points(DLIList<CubitPoint*> &point_list )
316 : : {
317 [ # # ]: 0 : for ( int i = 0; i < 4; i++ )
318 [ # # ]: 0 : point_list.append(this->point(i));
319 : 0 : }
320 : :
321 : :
322 : : //===========================================================================
323 : : // Function: facets
324 : : // Purpose: append the two facets comprising this quad. Virtual to FacetEntity
325 : : // Date: 8/1/03
326 : : // Author: sjowen
327 : : //===========================================================================
328 : 0 : void CubitQuadFacet::facets(DLIList<CubitFacet*> &facet_list )
329 : : {
330 [ # # ]: 0 : facet_list.append( this->get_tri_facet(0) );
331 [ # # ]: 0 : facet_list.append( this->get_tri_facet(1) );
332 : 0 : }
333 : :
334 : :
335 : : //===========================================================================
336 : : // Function: edges
337 : : // Purpose: append the edges of this quad. Virtual to FacetEntity
338 : : // Date: 8/1/03
339 : : // Author: sjowen
340 : : //===========================================================================
341 : 0 : void CubitQuadFacet::edges(DLIList<CubitFacetEdge*> &edge_list )
342 : : {
343 [ # # ]: 0 : for ( int i = 0; i < 4; i++ )
344 [ # # ]: 0 : edge_list.append(this->edge(i));
345 : 0 : }
346 : :
347 : : //===========================================================================
348 : : // Function: init_patch
349 : : // Purpose: computes the interior control points for the quad facet and
350 : : // stores them with the quad facet. Assumes edge control points
351 : : // on the four boundary edges have already been computed
352 : : // Date: 10/03
353 : : // Author: sjowen
354 : : //===========================================================================
355 : 0 : CubitStatus CubitQuadFacet::init_patch( )
356 : : {
357 : : // compute control point for the internal diagonal edge
358 : :
359 [ # # ][ # # ]: 0 : CubitFacetEdge *edge_ptr = this->point(0)->shared_edge( this->point(2) );
[ # # ]
360 [ # # ]: 0 : assert(edge_ptr != NULL);
361 [ # # ]: 0 : CubitPoint *pt0 = edge_ptr->point(0);
362 [ # # ]: 0 : CubitPoint *pt1 = edge_ptr->point(1);
363 [ # # ]: 0 : CubitVector P0 = pt0->coordinates();
364 [ # # ]: 0 : CubitVector P1 = pt1->coordinates();
365 [ # # ]: 0 : CubitVector N0 = pt0->normal( edge_ptr );
366 [ # # ]: 0 : CubitVector N1 = pt1->normal( edge_ptr );
367 [ # # ]: 0 : CubitVector T0 = P1 - P0;
368 [ # # ]: 0 : T0.normalize();
369 [ # # ]: 0 : CubitVector T1 = T0;
370 [ # # ][ # # ]: 0 : CubitVector Pi[3];
371 [ # # ]: 0 : CubitStatus stat = FacetEvalTool::init_edge_control_points( P0, P1, N0, N1, T0, T1, Pi);
372 [ # # ]: 0 : if (stat == CUBIT_FAILURE)
373 : 0 : return stat;
374 [ # # ]: 0 : edge_ptr->control_points( Pi, 4 );
375 : :
376 : : // compute control points on the triangles
377 : :
378 [ # # ][ # # ]: 0 : stat = this->get_tri_facet(0)->init_patch();
379 [ # # ]: 0 : if (stat != CUBIT_FAILURE)
380 [ # # ][ # # ]: 0 : stat = this->get_tri_facet(1)->init_patch();
381 : 0 : return stat;
382 : : }
383 : :
384 : :
385 : :
386 : : //EOF
|
