Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2006 Sandia National Laboratories. Developed at the
5 : : University of Wisconsin--Madison under SNL contract number
6 : : 624796. The U.S. Government and the University of Wisconsin
7 : : retain certain rights to this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : (2006) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file TriLagrangeShape.cpp
28 : : * \brief
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "TriLagrangeShape.hpp"
34 : : #include "MsqError.hpp"
35 : : #include <assert.h>
36 : :
37 : : namespace MBMesquite
38 : : {
39 : :
40 : 240850 : EntityTopology TriLagrangeShape::element_topology() const
41 : : {
42 : 240850 : return TRIANGLE;
43 : : }
44 : :
45 : 240849 : int TriLagrangeShape::num_nodes() const
46 : : {
47 : 240849 : return 6;
48 : : }
49 : :
50 : 1589658 : NodeSet TriLagrangeShape::sample_points( NodeSet ns ) const
51 : : {
52 [ + + ]: 1589658 : if( ns.have_any_mid_node() ) { ns.set_all_corner_nodes( TRIANGLE ); }
53 : : else
54 : : {
55 : 1579476 : ns.clear();
56 : 1579476 : ns.set_mid_face_node( 0 );
57 : : }
58 : 1589658 : return ns;
59 : : }
60 : :
61 : 3947 : void TriLagrangeShape::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out,
62 : : size_t& num_coeff, MsqError& err ) const
63 : : {
64 [ - + ]: 3947 : if( nodeset.have_any_mid_face_node() )
65 : : {
66 : : MSQ_SETERR( err )
67 [ # # ]: 0 : ( "TriLagrangeShape does not support mid-element nodes", MsqError::UNSUPPORTED_ELEMENT );
68 : 3947 : return;
69 : : }
70 : :
71 [ - + - - ]: 3947 : switch( loc.dimension )
72 : : {
73 : : case 0:
74 : 0 : num_coeff = 1;
75 : 0 : indices_out[0] = loc.number;
76 : 0 : coeff_out[0] = 1.0;
77 : 0 : break;
78 : : case 1:
79 [ - + ]: 3947 : if( nodeset.mid_edge_node( loc.number ) )
80 : : { // if mid-edge node is present
81 : 0 : num_coeff = 1;
82 : 0 : indices_out[0] = 3 + loc.number;
83 : 0 : coeff_out[0] = 1.0;
84 : : }
85 : : else
86 : : { // no mid node on edge
87 : 3947 : num_coeff = 2;
88 : 3947 : indices_out[0] = loc.number;
89 : 3947 : indices_out[1] = ( loc.number + 1 ) % 3;
90 : 3947 : coeff_out[0] = 0.5;
91 : 3947 : coeff_out[1] = 0.5;
92 : : }
93 : 3947 : break;
94 : : case 2:
95 : 0 : num_coeff = 3;
96 : 0 : indices_out[0] = 0;
97 : 0 : indices_out[1] = 1;
98 : 0 : indices_out[2] = 2;
99 : 0 : coeff_out[0] = 1.0 / 3.0;
100 : 0 : coeff_out[1] = 1.0 / 3.0;
101 : 0 : coeff_out[2] = 1.0 / 3.0;
102 [ # # ]: 0 : for( int i = 0; i < 3; ++i )
103 : : { // for each mid-edge node
104 [ # # ]: 0 : if( nodeset.mid_edge_node( i ) )
105 : : { // if node is present
106 : 0 : indices_out[num_coeff] = i + 3;
107 : : // coeff for mid-edge node
108 : 0 : coeff_out[num_coeff] = 4.0 / 9.0;
109 : : // adjust coeff for adj corner nodes
110 : 0 : coeff_out[i] -= 2.0 / 9.0;
111 : 0 : coeff_out[( i + 1 ) % 3] -= 2.0 / 9.0;
112 : : // update count
113 : 0 : ++num_coeff;
114 : : }
115 : : }
116 : 0 : break;
117 : : default:
118 : : MSQ_SETERR( err )
119 : : ( MsqError::UNSUPPORTED_ELEMENT,
120 : : "Request for dimension %d mapping function value"
121 : : "for a triangular element",
122 [ # # ]: 0 : loc.dimension );
123 : : }
124 : : }
125 : :
126 : 262633 : static inline void get_linear_derivatives( size_t* vertices, MsqVector< 2 >* derivs )
127 : : {
128 : 262633 : vertices[0] = 0;
129 : 262633 : vertices[1] = 1;
130 : 262633 : vertices[2] = 2;
131 : 262633 : derivs[0][0] = -1.0;
132 : 262633 : derivs[0][1] = -1.0;
133 : 262633 : derivs[1][0] = 1.0;
134 : 262633 : derivs[1][1] = 0.0;
135 : 262633 : derivs[2][0] = 0.0;
136 : 262633 : derivs[2][1] = 1.0;
137 : 262633 : }
138 : :
139 : 30546 : static void derivatives_at_corner( unsigned corner, NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs,
140 : : size_t& num_vtx )
141 : : {
142 : 30546 : num_vtx = 3;
143 : 30546 : get_linear_derivatives( vertices, derivs );
144 [ + + + - ]: 30546 : switch( corner )
145 : : {
146 : : case 0:
147 [ + - ]: 10182 : if( nodeset.mid_edge_node( 0 ) )
148 : : {
149 : 10182 : vertices[num_vtx] = 3;
150 : 10182 : derivs[num_vtx][0] = 4.0;
151 : 10182 : derivs[num_vtx][1] = 0.0;
152 : 10182 : ++num_vtx;
153 : 10182 : derivs[0][0] -= 2.0;
154 : 10182 : derivs[1][0] -= 2.0;
155 : : }
156 [ + + ]: 10182 : if( nodeset.mid_edge_node( 2 ) )
157 : : {
158 : 8454 : vertices[num_vtx] = 5;
159 : 8454 : derivs[num_vtx][0] = 0.0;
160 : 8454 : derivs[num_vtx][1] = 4.0;
161 : 8454 : ++num_vtx;
162 : 8454 : derivs[0][1] -= 2.0;
163 : 8454 : derivs[2][1] -= 2.0;
164 : : }
165 : 10182 : break;
166 : :
167 : : case 1:
168 [ + - ]: 10182 : if( nodeset.mid_edge_node( 0 ) )
169 : : {
170 : 10182 : vertices[num_vtx] = 3;
171 : 10182 : derivs[num_vtx][0] = -4.0;
172 : 10182 : derivs[num_vtx][1] = -4.0;
173 : 10182 : ++num_vtx;
174 : 10182 : derivs[0][0] += 2.0;
175 : 10182 : derivs[0][1] += 2.0;
176 : 10182 : derivs[1][0] += 2.0;
177 : 10182 : derivs[1][1] += 2.0;
178 : : }
179 [ + + ]: 10182 : if( nodeset.mid_edge_node( 1 ) )
180 : : {
181 : 8454 : vertices[num_vtx] = 4;
182 : 8454 : derivs[num_vtx][0] = 0.0;
183 : 8454 : derivs[num_vtx][1] = 4.0;
184 : 8454 : ++num_vtx;
185 : 8454 : derivs[1][1] -= 2.0;
186 : 8454 : derivs[2][1] -= 2.0;
187 : : }
188 : 10182 : break;
189 : :
190 : : case 2:
191 [ + + ]: 10182 : if( nodeset.mid_edge_node( 1 ) )
192 : : {
193 : 8454 : vertices[num_vtx] = 4;
194 : 8454 : derivs[num_vtx][0] = 4.0;
195 : 8454 : derivs[num_vtx][1] = 0.0;
196 : 8454 : ++num_vtx;
197 : 8454 : derivs[1][0] -= 2.0;
198 : 8454 : derivs[2][0] -= 2.0;
199 : : }
200 [ + + ]: 10182 : if( nodeset.mid_edge_node( 2 ) )
201 : : {
202 : 8454 : vertices[num_vtx] = 5;
203 : 8454 : derivs[num_vtx][0] = -4.0;
204 : 8454 : derivs[num_vtx][1] = -4.0;
205 : 8454 : ++num_vtx;
206 : 8454 : derivs[0][0] += 2.0;
207 : 8454 : derivs[0][1] += 2.0;
208 : 8454 : derivs[2][0] += 2.0;
209 : 8454 : derivs[2][1] += 2.0;
210 : : }
211 : 10182 : break;
212 : : }
213 : 30546 : }
214 : :
215 : : static const double edr[3][3] = { { 0.0, -2.0, 2.0 }, { 0.0, 2.0, 2.0 }, { 0.0, -2.0, -2.0 } };
216 : : static const double eds[3][3] = { { -2.0, -2.0, 0.0 }, { 2.0, 2.0, 0.0 }, { 2.0, -2.0, 0.0 } };
217 : :
218 : 27090 : static void derivatives_at_mid_edge( unsigned edge, NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs,
219 : : size_t& num_vtx )
220 : : {
221 : : // The mid-edge behavior is rather strange.
222 : : // A corner vertex contributes to the jacobian
223 : : // at the mid-point of the opposite edge unless
224 : : // one, but *not* both of the adjacent mid-edge
225 : : // nodes is present.
226 : :
227 : : // The value for each corner is incremented when
228 : : // a mid-side node is present. If the value is
229 : : // 0 when finished, the corner doesn't contribute.
230 : : // Initialize values to 0 for corners adjacent to
231 : : // edge so they are never zero.
232 : 27090 : int corner_count[3] = { 1, 1, 1 };
233 : 27090 : corner_count[( edge + 2 ) % 3] = -1;
234 : :
235 : : // begin with derivatives for linear tri
236 : 27090 : double corner_derivs[3][2] = { { -1.0, -1.0 }, { 1.0, 0.0 }, { 0.0, 1.0 } };
237 : :
238 : : // do mid-side nodes first
239 : 27090 : num_vtx = 0;
240 [ + + ]: 108360 : for( unsigned i = 0; i < 3; ++i )
241 : : {
242 [ + - ][ + + ]: 81270 : if( nodeset.mid_edge_node( i ) )
243 : : {
244 : 77814 : vertices[num_vtx] = i + 3;
245 [ + - ]: 77814 : derivs[num_vtx][0] = edr[i][edge];
246 [ + - ]: 77814 : derivs[num_vtx][1] = eds[i][edge];
247 : 77814 : ++num_vtx;
248 : :
249 : 77814 : int a = ( i + 1 ) % 3;
250 : 77814 : corner_derivs[i][0] -= 0.5 * edr[i][edge];
251 : 77814 : corner_derivs[i][1] -= 0.5 * eds[i][edge];
252 : 77814 : corner_derivs[a][0] -= 0.5 * edr[i][edge];
253 : 77814 : corner_derivs[a][1] -= 0.5 * eds[i][edge];
254 : 77814 : ++corner_count[i];
255 : 77814 : ++corner_count[a];
256 : : }
257 : : }
258 : :
259 : : // now add corner nodes to list
260 [ + + ]: 108360 : for( unsigned i = 0; i < 3; ++i )
261 : : {
262 [ + - ]: 81270 : if( corner_count[i] )
263 : : {
264 : 81270 : vertices[num_vtx] = i;
265 [ + - ]: 81270 : derivs[num_vtx][0] = corner_derivs[i][0];
266 [ + - ]: 81270 : derivs[num_vtx][1] = corner_derivs[i][1];
267 : 81270 : ++num_vtx;
268 : : }
269 : : }
270 : 27090 : }
271 : :
272 : : static const double fdr[] = { 0.0, 4.0 / 3.0, -4.0 / 3.0 };
273 : : static const double fds[] = { -4.0 / 3.0, 4.0 / 3.0, 0.0 };
274 : :
275 : 0 : static void derivatives_at_mid_elem( NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs, size_t& num_vtx )
276 : : {
277 : 0 : get_linear_derivatives( vertices, derivs );
278 : 0 : num_vtx = 3;
279 [ # # ]: 0 : for( unsigned i = 0; i < 3; ++i )
280 : : {
281 [ # # ]: 0 : if( nodeset.mid_edge_node( i ) )
282 : : {
283 : 0 : vertices[num_vtx] = i + 3;
284 : 0 : derivs[num_vtx][0] = fdr[i];
285 : 0 : derivs[num_vtx][1] = fds[i];
286 : 0 : ++num_vtx;
287 : :
288 : 0 : int a = ( i + 1 ) % 3;
289 : 0 : derivs[i][0] -= 0.5 * fdr[i];
290 : 0 : derivs[i][1] -= 0.5 * fds[i];
291 : 0 : derivs[a][0] -= 0.5 * fdr[i];
292 : 0 : derivs[a][1] -= 0.5 * fds[i];
293 : : }
294 : : }
295 : 0 : }
296 : :
297 : 289723 : void TriLagrangeShape::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out,
298 : : MsqVector< 2 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const
299 : : {
300 [ + + ]: 289723 : if( !nodeset.have_any_mid_node() )
301 : : {
302 : 232087 : num_vtx = 3;
303 : 232087 : get_linear_derivatives( vertex_indices_out, d_coeff_d_xi_out );
304 : 232087 : return;
305 : : }
306 : :
307 [ - + ]: 57636 : if( nodeset.have_any_mid_face_node() )
308 : : {
309 : : MSQ_SETERR( err )
310 [ # # ]: 0 : ( "TriLagrangeShape does not support mid-element nodes", MsqError::UNSUPPORTED_ELEMENT );
311 : 0 : return;
312 : : }
313 : :
314 [ + + - - ]: 57636 : switch( loc.dimension )
315 : : {
316 : : case 0:
317 : 30546 : derivatives_at_corner( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
318 : 30546 : break;
319 : : case 1:
320 : 27090 : derivatives_at_mid_edge( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
321 : 27090 : break;
322 : : case 2:
323 : 0 : derivatives_at_mid_elem( nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
324 : 0 : break;
325 : : default:
326 [ # # ]: 289723 : MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
327 : : }
328 : : }
329 : :
330 : 0 : void TriLagrangeShape::ideal( Sample, MsqMatrix< 3, 2 >& J, MsqError& ) const
331 : : {
332 : : // const double g = 1.074569931823542; // sqrt(2/sqrt(3));
333 : : // J(0,0) = g; J(0,1) = 0.5*g;
334 : : // J(1,0) = 0.0; J(1,1) = 1.0/g;
335 : : // J(2,0) = 0.0; J(2,1) = 0.0;
336 : 0 : const double a = 0.70710678118654746; // 1/sqrt(2)
337 : 0 : const double b = 1.3160740129524924; // 4th root of 3
338 : 0 : J( 0, 0 ) = -a / b;
339 : 0 : J( 0, 1 ) = a / b;
340 : 0 : J( 1, 0 ) = -a * b;
341 : 0 : J( 1, 1 ) = -a * b;
342 : 0 : J( 2, 0 ) = 0.0;
343 : 0 : J( 2, 1 ) = 0.0;
344 : 0 : }
345 : :
346 : : } // namespace MBMesquite
|
