Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2006 Lawrence Livermore National Laboratory. Under
5 : : the terms of Contract B545069 with the University of Wisconsin --
6 : : Madison, Lawrence Livermore National Laboratory retains certain
7 : : rights in 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 : : [email protected]
24 : :
25 : : ***************************************************************** */
26 : : /** \file QuadLagrangeShape.cpp
27 : : * \author Jason Kraftcheck
28 : : */
29 : :
30 : : #include "QuadLagrangeShape.hpp"
31 : : #include "LinearQuadrilateral.hpp"
32 : : #include "MsqError.hpp"
33 : :
34 : : namespace MBMesquite
35 : : {
36 : :
37 : 16501950 : EntityTopology QuadLagrangeShape::element_topology() const
38 : : {
39 : 16501950 : return QUADRILATERAL;
40 : : }
41 : :
42 : 13941450 : int QuadLagrangeShape::num_nodes() const
43 : : {
44 : 13941450 : return 9;
45 : : }
46 : :
47 : 708 : void QuadLagrangeShape::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out,
48 : : size_t& num_coeff, MsqError& err ) const
49 : : {
50 [ + - ]: 708 : if( !nodeset.have_any_mid_node() )
51 : : {
52 : 708 : LinearQuadrilateral::coefficients( loc, nodeset, coeff_out, indices_out, num_coeff );
53 : 708 : return;
54 : : }
55 : :
56 [ # # # # ]: 0 : switch( loc.dimension )
57 : : {
58 : : case 0:
59 : 0 : num_coeff = 1;
60 : 0 : indices_out[0] = loc.number;
61 : 0 : coeff_out[0] = 1.0;
62 : 0 : break;
63 : : case 1:
64 : 0 : coeff_out[0] = coeff_out[1] = coeff_out[2] = coeff_out[3] = coeff_out[4] = coeff_out[5] = coeff_out[6] =
65 : 0 : coeff_out[7] = coeff_out[8] = 0.0;
66 [ # # ]: 0 : if( nodeset.mid_edge_node( loc.number ) )
67 : : {
68 : : // if mid-edge node is present
69 : 0 : num_coeff = 1;
70 : 0 : indices_out[0] = loc.number + 4;
71 : 0 : coeff_out[0] = 1.0;
72 : : }
73 : : else
74 : : {
75 : : // If mid-edge node is not present, mapping function value
76 : : // for linear edge is even weight of adjacent vertices.
77 : 0 : num_coeff = 2;
78 : 0 : indices_out[0] = loc.number;
79 : 0 : indices_out[1] = ( loc.number + 1 ) % 4;
80 : 0 : coeff_out[0] = 0.5;
81 : 0 : coeff_out[1] = 0.5;
82 : : }
83 : 0 : break;
84 : : case 2:
85 [ # # ]: 0 : if( nodeset.mid_face_node( 0 ) )
86 : : { // if quad center node is present
87 : 0 : num_coeff = 1;
88 : 0 : indices_out[0] = 8;
89 : 0 : coeff_out[0] = 1.0;
90 : : }
91 : : else
92 : : {
93 : : // for linear element, (no mid-edge nodes), all corners contribute 1/4.
94 : 0 : num_coeff = 4;
95 : 0 : indices_out[0] = 0;
96 : 0 : indices_out[1] = 1;
97 : 0 : indices_out[2] = 2;
98 : 0 : indices_out[3] = 3;
99 : 0 : coeff_out[0] = 0.25;
100 : 0 : coeff_out[1] = 0.25;
101 : 0 : coeff_out[2] = 0.25;
102 : 0 : coeff_out[3] = 0.25;
103 : : // add in contribution for any mid-edge nodes present
104 [ # # ]: 0 : for( int i = 0; i < 4; ++i )
105 : : { // for each edge
106 [ # # ]: 0 : if( nodeset.mid_edge_node( i ) )
107 : : {
108 : 0 : indices_out[num_coeff] = i + 4;
109 : 0 : coeff_out[num_coeff] = 0.5;
110 : 0 : coeff_out[i] -= 0.25;
111 : 0 : coeff_out[( i + 1 ) % 4] -= 0.25;
112 : 0 : ++num_coeff;
113 : : }
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 quadrilateral element",
122 [ # # ]: 0 : loc.dimension );
123 : : }
124 : : }
125 : :
126 : 6625 : static void derivatives_at_corner( unsigned corner, NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs,
127 : : size_t& num_vtx )
128 : : {
129 : : static const unsigned xi_adj_corners[] = { 1, 0, 3, 2 };
130 : : static const unsigned xi_adj_edges[] = { 0, 0, 2, 2 };
131 : : static const unsigned eta_adj_corners[] = { 3, 2, 1, 0 };
132 : : static const unsigned eta_adj_edges[] = { 3, 1, 1, 3 };
133 : :
134 : : static const double corner_xi[] = { -1, 1, 1, -1 }; // xi values by corner
135 : : static const double corner_eta[] = { -1, -1, 1, 1 }; // eta values by corner
136 : : static const double other_xi[] = { 1, -1, -1, 1 }; // xi values for adjacent corner in xi direction
137 : : static const double other_eta[] = { 1, 1, -1, -1 }; // eta values for adjcent corner in eta direction
138 : : static const double mid_xi[] = { 2, -2, -2, 2 }; // xi values for mid-node in xi direction
139 : : static const double mid_eta[] = { 2, 2, -2, -2 }; // eta values for mid-node in eta direction
140 : :
141 : 6625 : num_vtx = 3;
142 : 6625 : vertices[0] = corner;
143 : 6625 : vertices[1] = xi_adj_corners[corner];
144 : 6625 : vertices[2] = eta_adj_corners[corner];
145 : :
146 : 6625 : derivs[0][0] = corner_xi[corner];
147 : 6625 : derivs[0][1] = corner_eta[corner];
148 : 6625 : derivs[1][0] = other_xi[corner];
149 : 6625 : derivs[1][1] = 0.0;
150 : 6625 : derivs[2][0] = 0.0;
151 : 6625 : derivs[2][1] = other_eta[corner];
152 : :
153 [ + - ]: 6625 : if( nodeset.mid_edge_node( xi_adj_edges[corner] ) )
154 : : {
155 : 6625 : vertices[num_vtx] = 4 + xi_adj_edges[corner];
156 : 6625 : derivs[num_vtx][0] = 2.0 * mid_xi[corner];
157 : 6625 : derivs[num_vtx][1] = 0.0;
158 : 6625 : derivs[0][0] -= mid_xi[corner];
159 : 6625 : derivs[1][0] -= mid_xi[corner];
160 : 6625 : ++num_vtx;
161 : : }
162 : :
163 [ + - ]: 6625 : if( nodeset.mid_edge_node( eta_adj_edges[corner] ) )
164 : : {
165 : 6625 : vertices[num_vtx] = 4 + eta_adj_edges[corner];
166 : 6625 : derivs[num_vtx][0] = 0.0;
167 : 6625 : derivs[num_vtx][1] = 2.0 * mid_eta[corner];
168 : 6625 : derivs[0][1] -= mid_eta[corner];
169 : 6625 : derivs[2][1] -= mid_eta[corner];
170 : 6625 : ++num_vtx;
171 : : }
172 : 6625 : }
173 : :
174 : 6596 : static void derivatives_at_mid_edge( unsigned edge, NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs,
175 : : size_t& num_vtx )
176 : : {
177 : : static const double values[][9] = { { -0.5, -0.5, 0.5, 0.5, -1.0, 2.0, 1.0, 2.0, 4.0 },
178 : : { -0.5, 0.5, 0.5, -0.5, -2.0, 1.0, -2.0, -1.0, -4.0 },
179 : : { -0.5, -0.5, 0.5, 0.5, -1.0, -2.0, 1.0, -2.0, -4.0 },
180 : : { -0.5, 0.5, 0.5, -0.5, 2.0, 1.0, 2.0, -1.0, 4.0 } };
181 : : static const double edge_values[][2] = { { -1, 1 }, { -1, 1 }, { 1, -1 }, { 1, -1 } };
182 : 6596 : const unsigned prev_corner = edge; // index of start vertex of edge
183 : 6596 : const unsigned next_corner = ( edge + 1 ) % 4; // index of end vertex of edge
184 : 6596 : const unsigned is_eta_edge = edge % 2; // edge is xi = +/- 0
185 : 6596 : const unsigned is_xi_edge = 1 - is_eta_edge; // edge is eta = +/- 0
186 : : // const unsigned mid_edge_node = edge + 4; // mid-edge node index
187 : 6596 : const unsigned prev_opposite = ( prev_corner + 3 ) % 4; // index of corner adjacent to prev_corner
188 : 6596 : const unsigned next_opposite = ( next_corner + 1 ) % 4; // index of corner adjacent to next_corner;
189 : :
190 : : // First do derivatives along edge (e.g. wrt xi if edge is eta = +/-1)
191 : 6596 : num_vtx = 2;
192 : 6596 : vertices[0] = prev_corner;
193 : 6596 : vertices[1] = next_corner;
194 [ + - ]: 6596 : derivs[0][is_eta_edge] = edge_values[edge][0];
195 [ + - ]: 6596 : derivs[0][is_xi_edge] = 0.0;
196 [ + - ]: 6596 : derivs[1][is_eta_edge] = edge_values[edge][1];
197 [ + - ]: 6596 : derivs[1][is_xi_edge] = 0.0;
198 : : // That's it for the edge-direction derivatives. No other vertices contribute.
199 : :
200 : : // Next handle the linear element case. Handle this as a special case first,
201 : : // so the generalized solution doesn't impact performance for linear elements
202 : : // too much.
203 [ + - ][ - + ]: 6596 : if( !nodeset.have_any_mid_node() )
204 : : {
205 : 0 : num_vtx = 4;
206 : 0 : vertices[2] = prev_opposite;
207 : 0 : vertices[3] = next_opposite;
208 [ # # ]: 0 : derivs[0][is_xi_edge] = values[edge][prev_corner];
209 [ # # ]: 0 : derivs[1][is_xi_edge] = values[edge][next_corner];
210 [ # # ]: 0 : derivs[2][is_xi_edge] = values[edge][prev_opposite];
211 [ # # ]: 0 : derivs[2][is_eta_edge] = 0.0;
212 [ # # ]: 0 : derivs[3][is_xi_edge] = values[edge][next_opposite];
213 [ # # ]: 0 : derivs[3][is_eta_edge] = 0.0;
214 : 6596 : return;
215 : : }
216 : :
217 : : // Initial (linear) contribution for each corner
218 : 6596 : double v[4] = { values[edge][0], values[edge][1], values[edge][2], values[edge][3] };
219 : :
220 : : // If mid-face node is present
221 : 6596 : double v8 = 0.0;
222 [ + - ][ - + ]: 6596 : if( nodeset.mid_face_node( 0 ) )
223 : : {
224 : 0 : v8 = values[edge][8];
225 : 0 : vertices[num_vtx] = 8;
226 [ # # ]: 0 : derivs[num_vtx][is_eta_edge] = 0.0;
227 [ # # ]: 0 : derivs[num_vtx][is_xi_edge] = v8;
228 : 0 : v[0] -= 0.25 * v8;
229 : 0 : v[1] -= 0.25 * v8;
230 : 0 : v[2] -= 0.25 * v8;
231 : 0 : v[3] -= 0.25 * v8;
232 : 0 : ++num_vtx;
233 : : }
234 : :
235 : : // If mid-edge nodes are present
236 [ + + ]: 32980 : for( unsigned i = 0; i < 4; ++i )
237 : : {
238 [ + - ][ + - ]: 26384 : if( nodeset.mid_edge_node( i ) )
239 : : {
240 : 26384 : const double value = values[edge][i + 4] - 0.5 * v8;
241 [ + - ]: 26384 : if( fabs( value ) > 0.125 )
242 : : {
243 : 26384 : v[i] -= 0.5 * value;
244 : 26384 : v[( i + 1 ) % 4] -= 0.5 * value;
245 : 26384 : vertices[num_vtx] = i + 4;
246 [ + - ]: 26384 : derivs[num_vtx][is_eta_edge] = 0.0;
247 [ + - ]: 26384 : derivs[num_vtx][is_xi_edge] = value;
248 : 26384 : ++num_vtx;
249 : : }
250 : : }
251 : : }
252 : :
253 : : // update values for adjacent corners
254 [ + - ]: 6596 : derivs[0][is_xi_edge] = v[prev_corner];
255 [ + - ]: 6596 : derivs[1][is_xi_edge] = v[next_corner];
256 : : // do other two corners
257 [ + - ]: 6596 : if( fabs( v[prev_opposite] ) > 0.125 )
258 : : {
259 : 6596 : vertices[num_vtx] = prev_opposite;
260 [ + - ]: 6596 : derivs[num_vtx][is_eta_edge] = 0.0;
261 [ + - ]: 6596 : derivs[num_vtx][is_xi_edge] = v[prev_opposite];
262 : 6596 : ++num_vtx;
263 : : }
264 [ + - ]: 6596 : if( fabs( v[next_opposite] ) > 0.125 )
265 : : {
266 : 6596 : vertices[num_vtx] = next_opposite;
267 [ + - ]: 6596 : derivs[num_vtx][is_eta_edge] = 0.0;
268 [ + - ]: 6596 : derivs[num_vtx][is_xi_edge] = v[next_opposite];
269 : 6596 : ++num_vtx;
270 : : }
271 : : }
272 : :
273 : 0 : static void derivatives_at_mid_elem( NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs, size_t& num_vtx )
274 : : {
275 : : // fast linear case
276 : : // This is provided as an optimization for linear elements.
277 : : // If this block of code were removed, the general-case code
278 : : // below should produce the same result.
279 [ # # ]: 0 : if( !nodeset.have_any_mid_node() )
280 : : {
281 : 0 : num_vtx = 4;
282 : 0 : vertices[0] = 0;
283 : 0 : derivs[0][0] = -0.5;
284 : 0 : derivs[0][1] = -0.5;
285 : 0 : vertices[1] = 1;
286 : 0 : derivs[1][0] = 0.5;
287 : 0 : derivs[1][1] = -0.5;
288 : 0 : vertices[2] = 2;
289 : 0 : derivs[2][0] = 0.5;
290 : 0 : derivs[2][1] = 0.5;
291 : 0 : vertices[3] = 3;
292 : 0 : derivs[3][0] = -0.5;
293 : 0 : derivs[3][1] = 0.5;
294 : 0 : return;
295 : : }
296 : :
297 : 0 : num_vtx = 0;
298 : :
299 : : // N_0
300 [ # # ]: 0 : if( !nodeset.both_edge_nodes( 0, 3 ) )
301 : : { // if eiter adjacent mid-edge node is missing
302 : 0 : vertices[num_vtx] = 0;
303 [ # # ]: 0 : derivs[num_vtx][0] = nodeset.mid_edge_node( 3 ) ? 0.0 : -0.5;
304 [ # # ]: 0 : derivs[num_vtx][1] = nodeset.mid_edge_node( 0 ) ? 0.0 : -0.5;
305 : 0 : ++num_vtx;
306 : : }
307 : :
308 : : // N_1
309 [ # # ]: 0 : if( !nodeset.both_edge_nodes( 0, 1 ) )
310 : : { // if eiter adjacent mid-edge node is missing
311 : 0 : vertices[num_vtx] = 1;
312 [ # # ]: 0 : derivs[num_vtx][0] = nodeset.mid_edge_node( 1 ) ? 0.0 : 0.5;
313 [ # # ]: 0 : derivs[num_vtx][1] = nodeset.mid_edge_node( 0 ) ? 0.0 : -0.5;
314 : 0 : ++num_vtx;
315 : : }
316 : :
317 : : // N_2
318 [ # # ]: 0 : if( !nodeset.both_edge_nodes( 1, 2 ) )
319 : : { // if eiter adjacent mid-edge node is missing
320 : 0 : vertices[num_vtx] = 2;
321 [ # # ]: 0 : derivs[num_vtx][0] = nodeset.mid_edge_node( 1 ) ? 0.0 : 0.5;
322 [ # # ]: 0 : derivs[num_vtx][1] = nodeset.mid_edge_node( 2 ) ? 0.0 : 0.5;
323 : 0 : ++num_vtx;
324 : : }
325 : :
326 : : // N_3
327 [ # # ]: 0 : if( !nodeset.both_edge_nodes( 2, 3 ) )
328 : : { // if eiter adjacent mid-edge node is missing
329 : 0 : vertices[num_vtx] = 3;
330 [ # # ]: 0 : derivs[num_vtx][0] = nodeset.mid_edge_node( 3 ) ? 0.0 : -0.5;
331 [ # # ]: 0 : derivs[num_vtx][1] = nodeset.mid_edge_node( 2 ) ? 0.0 : 0.5;
332 : 0 : ++num_vtx;
333 : : }
334 : :
335 : : // N_4
336 [ # # ]: 0 : if( nodeset.mid_edge_node( 0 ) )
337 : : {
338 : 0 : vertices[num_vtx] = 4;
339 : 0 : derivs[num_vtx][0] = 0.0;
340 : 0 : derivs[num_vtx][1] = -1.0;
341 : 0 : ++num_vtx;
342 : : }
343 : :
344 : : // N_5
345 [ # # ]: 0 : if( nodeset.mid_edge_node( 1 ) )
346 : : {
347 : 0 : vertices[num_vtx] = 5;
348 : 0 : derivs[num_vtx][0] = 1.0;
349 : 0 : derivs[num_vtx][1] = 0.0;
350 : 0 : ++num_vtx;
351 : : }
352 : :
353 : : // N_6
354 [ # # ]: 0 : if( nodeset.mid_edge_node( 2 ) )
355 : : {
356 : 0 : vertices[num_vtx] = 6;
357 : 0 : derivs[num_vtx][0] = 0.0;
358 : 0 : derivs[num_vtx][1] = 1.0;
359 : 0 : ++num_vtx;
360 : : }
361 : :
362 : : // N_7
363 [ # # ]: 0 : if( nodeset.mid_edge_node( 3 ) )
364 : : {
365 : 0 : vertices[num_vtx] = 7;
366 : 0 : derivs[num_vtx][0] = -1.0;
367 : 0 : derivs[num_vtx][1] = 0.0;
368 : 0 : ++num_vtx;
369 : : }
370 : :
371 : : // N_8 (mid-quad node) never contributes to Jacobian at element center!!!
372 : : }
373 : :
374 : 19387043 : void QuadLagrangeShape::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out,
375 : : MsqVector< 2 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const
376 : : {
377 [ + + ]: 19387043 : if( !nodeset.have_any_mid_node() )
378 : : {
379 : 19373822 : LinearQuadrilateral::derivatives( loc, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
380 : 19387043 : return;
381 : : }
382 : :
383 [ + + - - ]: 13221 : switch( loc.dimension )
384 : : {
385 : : case 0:
386 : 6625 : derivatives_at_corner( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
387 : 6625 : break;
388 : : case 1:
389 : 6596 : derivatives_at_mid_edge( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
390 : 6596 : break;
391 : : case 2:
392 : 0 : derivatives_at_mid_elem( nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
393 : 0 : break;
394 : : default:
395 [ # # ]: 0 : MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
396 : : }
397 : : }
398 : :
399 : 0 : void QuadLagrangeShape::ideal( Sample, MsqMatrix< 3, 2 >& J, MsqError& ) const
400 : : {
401 : 0 : J( 0, 0 ) = J( 1, 1 ) = 1.0;
402 : 0 : J( 0, 1 ) = J( 1, 0 ) = 0.0;
403 : 0 : J( 2, 0 ) = J( 2, 1 ) = 0.0;
404 : 0 : }
405 : :
406 : : } // namespace MBMesquite
|
