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 : : (2006) [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : /** \file LinearPrism.cpp
28 : : * \brief mapping function for linear prism
29 : : * \author Jason Kraftcheck
30 : : */
31 : :
32 : : #include "Mesquite.hpp"
33 : : #include "MsqError.hpp"
34 : : #include "LinearPrism.hpp"
35 : :
36 : : namespace MBMesquite
37 : : {
38 : :
39 : : static const char* nonlinear_error = "Attempt to use LinearPrism mapping function for a nonlinear element\n";
40 : :
41 : 0 : EntityTopology LinearPrism::element_topology() const
42 : : {
43 : 0 : return PRISM;
44 : : }
45 : :
46 : 0 : int LinearPrism::num_nodes() const
47 : : {
48 : 0 : return 6;
49 : : }
50 : :
51 : : static const int edge_beg[] = { 0, 1, 2, 0, 1, 2, 3, 4, 5 };
52 : : static const int edge_end[] = { 1, 2, 0, 3, 4, 5, 4, 5, 3 };
53 : : static const int faces[5][5] = { { 4, 0, 1, 4, 3 },
54 : : { 4, 1, 2, 5, 4 },
55 : : { 4, 2, 0, 3, 5 },
56 : : { 3, 0, 1, 2, -1 },
57 : : { 3, 3, 4, 5, -1 } };
58 : :
59 : 0 : static void coefficients_at_corner( unsigned corner, double* coeff_out, size_t* indices_out, size_t& num_coeff )
60 : : {
61 : 0 : num_coeff = 1;
62 : 0 : indices_out[0] = corner;
63 : 0 : coeff_out[0] = 1.0;
64 : 0 : }
65 : :
66 : 0 : static void coefficients_at_mid_edge( unsigned edge, double* coeff_out, size_t* indices_out, size_t& num_coeff )
67 : : {
68 : 0 : num_coeff = 2;
69 : 0 : indices_out[0] = edge_beg[edge];
70 : 0 : indices_out[1] = edge_end[edge];
71 : 0 : coeff_out[0] = 0.5;
72 : 0 : coeff_out[1] = 0.5;
73 : 0 : }
74 : :
75 : 0 : static void coefficients_at_mid_face( unsigned face, double* coeff_out, size_t* indices_out, size_t& num_coeff )
76 : : {
77 : : double f;
78 [ # # ]: 0 : if( faces[face][0] == 4 )
79 : : {
80 : 0 : num_coeff = 4;
81 : 0 : f = 0.25;
82 : 0 : indices_out[3] = faces[face][4];
83 : 0 : coeff_out[3] = f;
84 : : }
85 : : else
86 : : {
87 : 0 : num_coeff = 3;
88 : 0 : f = MSQ_ONE_THIRD;
89 : : }
90 : :
91 : 0 : coeff_out[0] = f;
92 : 0 : coeff_out[1] = f;
93 : 0 : coeff_out[2] = f;
94 : 0 : indices_out[0] = faces[face][1];
95 : 0 : indices_out[1] = faces[face][2];
96 : 0 : indices_out[2] = faces[face][3];
97 : 0 : }
98 : :
99 : 0 : static void coefficients_at_mid_elem( double* coeff_out, size_t* indices_out, size_t& num_coeff )
100 : : {
101 : 0 : num_coeff = 6;
102 : 0 : const double sixth = 1.0 / 6.0;
103 : 0 : coeff_out[0] = sixth;
104 : 0 : coeff_out[1] = sixth;
105 : 0 : coeff_out[2] = sixth;
106 : 0 : coeff_out[3] = sixth;
107 : 0 : coeff_out[4] = sixth;
108 : 0 : coeff_out[5] = sixth;
109 : 0 : indices_out[0] = 0;
110 : 0 : indices_out[1] = 1;
111 : 0 : indices_out[2] = 2;
112 : 0 : indices_out[3] = 3;
113 : 0 : indices_out[4] = 4;
114 : 0 : indices_out[5] = 5;
115 : 0 : }
116 : :
117 : 0 : void LinearPrism::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out, size_t& num_coeff,
118 : : MsqError& err ) const
119 : : {
120 [ # # ]: 0 : if( nodeset.have_any_mid_node() )
121 : : {
122 [ # # ]: 0 : MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
123 : 0 : return;
124 : : }
125 : :
126 [ # # # # : 0 : switch( loc.dimension )
# ]
127 : : {
128 : : case 0:
129 : 0 : coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
130 : 0 : break;
131 : : case 1:
132 : 0 : coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff );
133 : 0 : break;
134 : : case 2:
135 : 0 : coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff );
136 : 0 : break;
137 : : case 3:
138 : 0 : coefficients_at_mid_elem( coeff_out, indices_out, num_coeff );
139 : 0 : break;
140 : : default:
141 [ # # ]: 0 : MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
142 : : }
143 : : }
144 : :
145 : 0 : static void derivatives_at_corner( unsigned corner, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out,
146 : : size_t& num_vtx )
147 : : {
148 : 0 : int tri = ( corner / 3 ); // 0 for xi=0, 1 for xi=1
149 : 0 : int tv = corner % 3; // index of corner with xi=constant triangle
150 : :
151 : 0 : num_vtx = 4;
152 : : // three vertices within the xi=constant triangle
153 : 0 : vertex_indices_out[0] = 3 * tri;
154 : 0 : vertex_indices_out[1] = 3 * tri + 1;
155 : 0 : vertex_indices_out[2] = 3 * tri + 2;
156 : : // vertex adjacent to corner in other triangle
157 : 0 : vertex_indices_out[3] = 3 - 6 * tri + corner;
158 : :
159 : : // three vertices within the xi=constant triangle
160 : 0 : d_coeff_d_xi_out[0][0] = 0.0;
161 : 0 : d_coeff_d_xi_out[0][1] = -1.0;
162 : 0 : d_coeff_d_xi_out[0][2] = -1.0;
163 : 0 : d_coeff_d_xi_out[1][0] = 0.0;
164 : 0 : d_coeff_d_xi_out[1][1] = 1.0;
165 : 0 : d_coeff_d_xi_out[1][2] = 0.0;
166 : 0 : d_coeff_d_xi_out[2][0] = 0.0;
167 : 0 : d_coeff_d_xi_out[2][1] = 0.0;
168 : 0 : d_coeff_d_xi_out[2][2] = 1.0;
169 : : // fix dxi value for input corner
170 : 0 : d_coeff_d_xi_out[tv][0] = 2 * tri - 1;
171 : : // vertex adjacent to corner in other triangle
172 : 0 : d_coeff_d_xi_out[3][0] = 1 - 2 * tri;
173 : 0 : d_coeff_d_xi_out[3][1] = 0.0;
174 : 0 : d_coeff_d_xi_out[3][2] = 0.0;
175 : 0 : }
176 : :
177 : 0 : static void derivatives_at_mid_edge( unsigned edge, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out,
178 : : size_t& num_vtx )
179 : : {
180 : : int opp; // vertex opposite edge in same triagle
181 : :
182 [ # # # # ]: 0 : switch( edge / 3 )
183 : : {
184 : : case 0: // triangle at xi = 0
185 : 0 : opp = ( edge + 2 ) % 3;
186 : :
187 : 0 : num_vtx = 5;
188 : : // vertices in this xi = 0 triagnle
189 : 0 : vertex_indices_out[0] = 0;
190 : 0 : vertex_indices_out[1] = 1;
191 : 0 : vertex_indices_out[2] = 2;
192 : : // adjacent vertices in xi = 1 triangle
193 : 0 : vertex_indices_out[3] = 3 + edge;
194 : 0 : vertex_indices_out[4] = 3 + ( edge + 1 ) % 3;
195 : :
196 : : // vertices in this xi = 0 triagnle
197 : 0 : d_coeff_d_xi_out[0][0] = -0.5;
198 : 0 : d_coeff_d_xi_out[0][1] = -1.0;
199 : 0 : d_coeff_d_xi_out[0][2] = -1.0;
200 : 0 : d_coeff_d_xi_out[1][0] = -0.5;
201 : 0 : d_coeff_d_xi_out[1][1] = 1.0;
202 : 0 : d_coeff_d_xi_out[1][2] = 0.0;
203 : 0 : d_coeff_d_xi_out[2][0] = -0.5;
204 : 0 : d_coeff_d_xi_out[2][1] = 0.0;
205 : 0 : d_coeff_d_xi_out[2][2] = 1.0;
206 : : // clear dxi for vertex opposite edge in xi = 0 triangle
207 : 0 : d_coeff_d_xi_out[opp][0] = 0.0;
208 : : // adjacent vertices in xi = 1 triangle
209 : 0 : d_coeff_d_xi_out[3][0] = 0.5;
210 : 0 : d_coeff_d_xi_out[3][1] = 0.0;
211 : 0 : d_coeff_d_xi_out[3][2] = 0.0;
212 : 0 : d_coeff_d_xi_out[4][0] = 0.5;
213 : 0 : d_coeff_d_xi_out[4][1] = 0.0;
214 : 0 : d_coeff_d_xi_out[4][2] = 0.0;
215 : 0 : break;
216 : :
217 : : case 1: // lateral edges (not in either triangle)
218 : 0 : num_vtx = 6;
219 : 0 : vertex_indices_out[0] = 0;
220 : 0 : vertex_indices_out[1] = 1;
221 : 0 : vertex_indices_out[2] = 2;
222 : 0 : vertex_indices_out[3] = 3;
223 : 0 : vertex_indices_out[4] = 4;
224 : 0 : vertex_indices_out[5] = 5;
225 : :
226 : : // set all deta & dzeta values, zero all dxi values
227 : 0 : d_coeff_d_xi_out[0][0] = 0.0;
228 : 0 : d_coeff_d_xi_out[0][1] = -0.5;
229 : 0 : d_coeff_d_xi_out[0][2] = -0.5;
230 : 0 : d_coeff_d_xi_out[1][0] = 0.0;
231 : 0 : d_coeff_d_xi_out[1][1] = 0.5;
232 : 0 : d_coeff_d_xi_out[1][2] = 0.0;
233 : 0 : d_coeff_d_xi_out[2][0] = 0.0;
234 : 0 : d_coeff_d_xi_out[2][1] = 0.0;
235 : 0 : d_coeff_d_xi_out[2][2] = 0.5;
236 : 0 : d_coeff_d_xi_out[3][0] = 0.0;
237 : 0 : d_coeff_d_xi_out[3][1] = -0.5;
238 : 0 : d_coeff_d_xi_out[3][2] = -0.5;
239 : 0 : d_coeff_d_xi_out[4][0] = 0.0;
240 : 0 : d_coeff_d_xi_out[4][1] = 0.5;
241 : 0 : d_coeff_d_xi_out[4][2] = 0.0;
242 : 0 : d_coeff_d_xi_out[5][0] = 0.0;
243 : 0 : d_coeff_d_xi_out[5][1] = 0.0;
244 : 0 : d_coeff_d_xi_out[5][2] = 0.5;
245 : :
246 : : // set dxi values for end points of edge
247 : 0 : d_coeff_d_xi_out[( edge - 3 )][0] = -1;
248 : 0 : d_coeff_d_xi_out[edge][0] = 1;
249 : 0 : break;
250 : :
251 : : case 2: // triangle at xi = 1
252 : 0 : opp = ( edge + 2 ) % 3;
253 : :
254 : 0 : num_vtx = 5;
255 : : // vertices in this xi = 1 triagnle
256 : 0 : vertex_indices_out[0] = 3;
257 : 0 : vertex_indices_out[1] = 4;
258 : 0 : vertex_indices_out[2] = 5;
259 : : // adjacent vertices in xi = 1 triangle
260 : 0 : vertex_indices_out[3] = edge - 6;
261 : 0 : vertex_indices_out[4] = ( edge - 5 ) % 3;
262 : :
263 : : // vertices in this xi = 1 triagnle
264 : 0 : d_coeff_d_xi_out[0][0] = 0.5;
265 : 0 : d_coeff_d_xi_out[0][1] = -1.0;
266 : 0 : d_coeff_d_xi_out[0][2] = -1.0;
267 : 0 : d_coeff_d_xi_out[1][0] = 0.5;
268 : 0 : d_coeff_d_xi_out[1][1] = 1.0;
269 : 0 : d_coeff_d_xi_out[1][2] = 0.0;
270 : 0 : d_coeff_d_xi_out[2][0] = 0.5;
271 : 0 : d_coeff_d_xi_out[2][1] = 0.0;
272 : 0 : d_coeff_d_xi_out[2][2] = 1.0;
273 : : // clear dxi for vertex opposite edge in xi = 1 triangle
274 : 0 : d_coeff_d_xi_out[opp][0] = 0.0;
275 : : // adjacent vertices in xi = 0 triangle
276 : 0 : d_coeff_d_xi_out[3][0] = -0.5;
277 : 0 : d_coeff_d_xi_out[3][1] = 0.0;
278 : 0 : d_coeff_d_xi_out[3][2] = 0.0;
279 : 0 : d_coeff_d_xi_out[4][0] = -0.5;
280 : 0 : d_coeff_d_xi_out[4][1] = 0.0;
281 : 0 : d_coeff_d_xi_out[4][2] = 0.0;
282 : 0 : break;
283 : : }
284 : 0 : }
285 : 0 : static void derivatives_at_mid_face( unsigned face, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out,
286 : : size_t& num_vtx )
287 : : {
288 : 0 : num_vtx = 6;
289 : 0 : vertex_indices_out[0] = 0;
290 : 0 : vertex_indices_out[1] = 1;
291 : 0 : vertex_indices_out[2] = 2;
292 : 0 : vertex_indices_out[3] = 3;
293 : 0 : vertex_indices_out[4] = 4;
294 : 0 : vertex_indices_out[5] = 5;
295 : :
296 : : int opp; // start vtx of edge opposite from quad face
297 : : int tri_offset; // offset in d_coeff_d_xi_out for triangle containing edge
298 : :
299 [ # # ]: 0 : if( face < 3 )
300 : : { // quad face
301 : : // set all values
302 : 0 : d_coeff_d_xi_out[0][0] = -0.5;
303 : 0 : d_coeff_d_xi_out[0][1] = -0.5;
304 : 0 : d_coeff_d_xi_out[0][2] = -0.5;
305 : 0 : d_coeff_d_xi_out[1][0] = -0.5;
306 : 0 : d_coeff_d_xi_out[1][1] = 0.5;
307 : 0 : d_coeff_d_xi_out[1][2] = 0.0;
308 : 0 : d_coeff_d_xi_out[2][0] = -0.5;
309 : 0 : d_coeff_d_xi_out[2][1] = 0.0;
310 : 0 : d_coeff_d_xi_out[2][2] = 0.5;
311 : 0 : d_coeff_d_xi_out[3][0] = 0.5;
312 : 0 : d_coeff_d_xi_out[3][1] = -0.5;
313 : 0 : d_coeff_d_xi_out[3][2] = -0.5;
314 : 0 : d_coeff_d_xi_out[4][0] = 0.5;
315 : 0 : d_coeff_d_xi_out[4][1] = 0.5;
316 : 0 : d_coeff_d_xi_out[4][2] = 0.0;
317 : 0 : d_coeff_d_xi_out[5][0] = 0.5;
318 : 0 : d_coeff_d_xi_out[5][1] = 0.0;
319 : 0 : d_coeff_d_xi_out[5][2] = 0.5;
320 : : // clear dxi for ends of edge opposite from face
321 : 0 : opp = ( face + 2 ) % 3;
322 : 0 : d_coeff_d_xi_out[opp][0] = 0.0;
323 : 0 : d_coeff_d_xi_out[( opp + 3 )][0] = 0.0;
324 : : }
325 : : else
326 : : { // triangular faces
327 : : // set all xi values, zero all other values
328 : 0 : const double third = 1. / 3;
329 : 0 : d_coeff_d_xi_out[0][0] = -third;
330 : 0 : d_coeff_d_xi_out[0][1] = 0;
331 : 0 : d_coeff_d_xi_out[0][2] = 0;
332 : 0 : d_coeff_d_xi_out[1][0] = -third;
333 : 0 : d_coeff_d_xi_out[1][1] = 0;
334 : 0 : d_coeff_d_xi_out[1][2] = 0;
335 : 0 : d_coeff_d_xi_out[2][0] = -third;
336 : 0 : d_coeff_d_xi_out[2][1] = 0;
337 : 0 : d_coeff_d_xi_out[2][2] = 0;
338 : 0 : d_coeff_d_xi_out[3][0] = third;
339 : 0 : d_coeff_d_xi_out[3][1] = 0;
340 : 0 : d_coeff_d_xi_out[3][2] = 0;
341 : 0 : d_coeff_d_xi_out[4][0] = third;
342 : 0 : d_coeff_d_xi_out[4][1] = 0;
343 : 0 : d_coeff_d_xi_out[4][2] = 0;
344 : 0 : d_coeff_d_xi_out[5][0] = third;
345 : 0 : d_coeff_d_xi_out[5][1] = 0;
346 : 0 : d_coeff_d_xi_out[5][2] = 0;
347 : : // set deta and dzeta values for vertices in same triangle as edge
348 : 0 : tri_offset = 3 * ( face - 3 ); // either 0 or 3
349 : 0 : d_coeff_d_xi_out[tri_offset][1] = -1.0;
350 : 0 : d_coeff_d_xi_out[tri_offset][2] = -1.0;
351 : 0 : d_coeff_d_xi_out[tri_offset + 1][1] = 1.0;
352 : 0 : d_coeff_d_xi_out[tri_offset + 2][2] = 1.0;
353 : : }
354 : 0 : }
355 : 0 : static void derivatives_at_mid_elem( size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx )
356 : : {
357 : 0 : const double third = 1. / 3;
358 : :
359 : 0 : num_vtx = 6;
360 : : ;
361 : 0 : vertex_indices_out[0] = 0;
362 : 0 : vertex_indices_out[1] = 1;
363 : 0 : vertex_indices_out[2] = 2;
364 : 0 : vertex_indices_out[3] = 3;
365 : 0 : vertex_indices_out[4] = 4;
366 : 0 : vertex_indices_out[5] = 5;
367 : :
368 : 0 : d_coeff_d_xi_out[0][0] = -third;
369 : 0 : d_coeff_d_xi_out[0][1] = -0.5;
370 : 0 : d_coeff_d_xi_out[0][2] = -0.5;
371 : 0 : d_coeff_d_xi_out[1][0] = -third;
372 : 0 : d_coeff_d_xi_out[1][1] = 0.5;
373 : 0 : d_coeff_d_xi_out[1][2] = 0.0;
374 : 0 : d_coeff_d_xi_out[2][0] = -third;
375 : 0 : d_coeff_d_xi_out[2][1] = 0.0;
376 : 0 : d_coeff_d_xi_out[2][2] = 0.5;
377 : 0 : d_coeff_d_xi_out[3][0] = third;
378 : 0 : d_coeff_d_xi_out[3][1] = -0.5;
379 : 0 : d_coeff_d_xi_out[3][2] = -0.5;
380 : 0 : d_coeff_d_xi_out[4][0] = third;
381 : 0 : d_coeff_d_xi_out[4][1] = 0.5;
382 : 0 : d_coeff_d_xi_out[4][2] = 0.0;
383 : 0 : d_coeff_d_xi_out[5][0] = third;
384 : 0 : d_coeff_d_xi_out[5][1] = 0.0;
385 : 0 : d_coeff_d_xi_out[5][2] = 0.5;
386 : 0 : }
387 : :
388 : 0 : void LinearPrism::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out,
389 : : MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const
390 : : {
391 [ # # ]: 0 : if( nodeset.have_any_mid_node() )
392 : : {
393 [ # # ]: 0 : MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
394 : 0 : return;
395 : : }
396 : :
397 [ # # # # : 0 : switch( loc.dimension )
# ]
398 : : {
399 : : case 0:
400 : 0 : derivatives_at_corner( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
401 : 0 : break;
402 : : case 1:
403 : 0 : derivatives_at_mid_edge( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
404 : 0 : break;
405 : : case 2:
406 : 0 : derivatives_at_mid_face( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
407 : 0 : break;
408 : : case 3:
409 : 0 : derivatives_at_mid_elem( vertex_indices_out, d_coeff_d_xi_out, num_vtx );
410 : 0 : break;
411 : : default:
412 [ # # ]: 0 : MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
413 : : }
414 : : }
415 : :
416 : 0 : void LinearPrism::ideal( Sample, MsqMatrix< 3, 3 >& J, MsqError& ) const
417 : : {
418 : 0 : const double a = 0.52455753171082409; // 2^(-2/3) * 3^(-1/6)
419 : 0 : const double b = 0.90856029641606983; // a * sqrt(3) = 1/2 cbrt(6)
420 : :
421 : 0 : J( 0, 0 ) = 2 * a;
422 : 0 : J( 0, 1 ) = 0.0;
423 : 0 : J( 0, 2 ) = 0.0;
424 : 0 : J( 1, 0 ) = 0.0;
425 : 0 : J( 1, 1 ) = 2 * a;
426 : 0 : J( 1, 2 ) = a;
427 : 0 : J( 2, 0 ) = 0.0;
428 : 0 : J( 2, 1 ) = 0.0;
429 : 0 : J( 2, 2 ) = b;
430 : 0 : }
431 : :
432 : : } // namespace MBMesquite
|
