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