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 LinearPyramid.cpp
29 * \brief LinearPyramid implementation
30 * \author Jason Kraftcheck
31 */
32
33 #include "Mesquite.hpp"
34 #include "LinearPyramid.hpp"
35 #include "MsqError.hpp"
36
37 namespace MBMesquite {
38
39 static const char* nonlinear_error
40 = "Attempt to use LinearTriangle mapping function for a nonlinear element\n";
41
set_equal_derivatives(double value,size_t * indices,MsqVector<3> * derivs,size_t & num_vtx)42 static inline void set_equal_derivatives( double value,
43 size_t* indices,
44 MsqVector<3>* derivs,
45 size_t& num_vtx )
46 {
47 num_vtx = 5;
48 indices[0] = 0;
49 indices[1] = 1;
50 indices[2] = 2;
51 indices[3] = 3;
52 indices[4] = 4;
53
54 derivs[0][0] = -value;
55 derivs[0][1] = -value;
56 derivs[0][2] = -0.25;
57
58 derivs[1][0] = value;
59 derivs[1][1] = -value;
60 derivs[1][2] = -0.25;
61
62 derivs[2][0] = value;
63 derivs[2][1] = value;
64 derivs[2][2] = -0.25;
65
66 derivs[3][0] = -value;
67 derivs[3][1] = value;
68 derivs[3][2] = -0.25;
69
70 derivs[4][0] = 0.0;
71 derivs[4][1] = 0.0;
72 derivs[4][2] = 1.0;
73 }
74
set_edge_derivatives(unsigned base_corner,double value,size_t * indices,MsqVector<3> * derivs,size_t & num_vtx)75 static inline void set_edge_derivatives( unsigned base_corner,
76 double value,
77 size_t* indices,
78 MsqVector<3>* derivs,
79 size_t& num_vtx )
80 {
81 const int direction = base_corner % 2;
82 const int edge_beg = base_corner;
83 const int edge_end = (base_corner+1)%4;
84 const int adj_end = (base_corner+2)%4;
85 const int adj_beg = (base_corner+3)%4;
86 const int dir_sign = 2*(edge_beg/2) - 1;
87 const int oth_sign = 2*((edge_beg+1)/2%2) - 1;
88
89 num_vtx = 5;
90 indices[0] = edge_beg;
91 indices[1] = edge_end;
92 indices[2] = adj_end;
93 indices[3] = adj_beg;
94 indices[4] = 4;
95
96 derivs[0][ direction] = 2 * dir_sign * value;
97 derivs[0][1-direction] = oth_sign * value;
98 derivs[0][2] = -0.5;
99
100 derivs[1][ direction] = -2 * dir_sign * value;
101 derivs[1][1-direction] = oth_sign * value;
102 derivs[1][2] = -0.5;
103
104 derivs[2][ direction] = 0.0;
105 derivs[2][1-direction] = -oth_sign * value;
106 derivs[2][2] = 0.0;
107
108 derivs[3][ direction] = 0.0;
109 derivs[3][1-direction] = -oth_sign * value;
110 derivs[3][2] = 0.0;
111
112 derivs[4][0] = 0.0;
113 derivs[4][1] = 0.0;
114 derivs[4][2] = 1.0;
115 }
116
set_corner_derivatives(unsigned corner,double value,size_t * indices,MsqVector<3> * derivs,size_t & num_vtx)117 static inline void set_corner_derivatives( unsigned corner,
118 double value,
119 size_t* indices,
120 MsqVector<3>* derivs,
121 size_t& num_vtx )
122 {
123 const unsigned adj_in_xi = (5 - corner) % 4;
124 const unsigned adj_in_eta = 3 - corner;
125
126 const int dxi_sign = 2*((corner+1)/2%2)-1;
127 const int deta_sign = 2*(corner/2) - 1;
128 const double dxi_value = dxi_sign * value;
129 const double deta_value = deta_sign * value;
130
131 num_vtx = 4;
132 indices[0] = corner;
133 indices[1] = adj_in_xi;
134 indices[2] = adj_in_eta;
135 indices[3] = 4;
136
137 derivs[0][0] = dxi_value;
138 derivs[0][1] = deta_value;
139 derivs[0][2] = -1.0;
140
141 derivs[1][0] = -dxi_value;
142 derivs[1][1] = 0.0;
143 derivs[1][2] = 0.0;
144
145 derivs[2][0] = 0.0;
146 derivs[2][1] = -deta_value;
147 derivs[2][2] = 0.0;
148
149 derivs[3][0] = 0.0;
150 derivs[3][1] = 0.0;
151 derivs[3][2] = 1.0;
152 }
153
element_topology() const154 EntityTopology LinearPyramid::element_topology() const
155 { return PYRAMID; }
156
num_nodes() const157 int LinearPyramid::num_nodes() const
158 { return 5; }
159
sample_points(NodeSet) const160 NodeSet LinearPyramid::sample_points( NodeSet ) const
161 {
162 NodeSet result;
163 result.set_all_corner_nodes(PYRAMID);
164 result.clear_corner_node(4);
165 return result;
166 }
167
coefficients_at_corner(unsigned corner,double * coeff_out,size_t * indices_out,size_t & num_coeff)168 static void coefficients_at_corner( unsigned corner,
169 double* coeff_out,
170 size_t* indices_out,
171 size_t& num_coeff )
172 {
173 num_coeff = 1;
174 indices_out[0] = corner;
175 coeff_out[0] = 1.0;
176 }
177
178
coefficients_at_mid_edge(unsigned edge,double * coeff_out,size_t * indices_out,size_t & num_coeff)179 static void coefficients_at_mid_edge( unsigned edge,
180 double* coeff_out,
181 size_t* indices_out,
182 size_t& num_coeff )
183 {
184 num_coeff = 2;
185 coeff_out[0] = 0.5;
186 coeff_out[1] = 0.5;
187
188 if (edge < 4) {
189 indices_out[0] = edge;
190 indices_out[1] = (edge+1)%4;
191 }
192 else {
193 indices_out[0] = edge-4;
194 indices_out[1] = 4;
195 }
196 }
197
coefficients_at_mid_face(unsigned face,double * coeff_out,size_t * indices_out,size_t & num_coeff)198 static void coefficients_at_mid_face( unsigned face,
199 double* coeff_out,
200 size_t* indices_out,
201 size_t& num_coeff )
202 {
203 if (face == 4) {
204 num_coeff = 4;
205 coeff_out[0] = 0.25;
206 coeff_out[1] = 0.25;
207 coeff_out[2] = 0.25;
208 coeff_out[3] = 0.25;
209 indices_out[0] = 0;
210 indices_out[1] = 1;
211 indices_out[2] = 2;
212 indices_out[3] = 3;
213 }
214 else {
215 num_coeff = 3;
216 indices_out[0] = face;
217 indices_out[1] = (face+1)%4;
218 indices_out[2] = 4;
219 coeff_out[0] = 0.25;
220 coeff_out[1] = 0.25;
221 coeff_out[2] = 0.50;
222 }
223 }
224
coefficients_at_mid_elem(double * coeff_out,size_t * indices_out,size_t & num_coeff)225 static void coefficients_at_mid_elem( double* coeff_out,
226 size_t* indices_out,
227 size_t& num_coeff )
228 {
229 num_coeff = 5;
230 coeff_out[0] = 0.125;
231 coeff_out[1] = 0.125;
232 coeff_out[2] = 0.125;
233 coeff_out[3] = 0.125;
234 coeff_out[4] = 0.500;
235 indices_out[0] = 0;
236 indices_out[1] = 1;
237 indices_out[2] = 2;
238 indices_out[3] = 3;
239 indices_out[4] = 4;
240 }
241
coefficients(Sample loc,NodeSet nodeset,double * coeff_out,size_t * indices_out,size_t & num_coeff,MsqError & err) const242 void LinearPyramid::coefficients( Sample loc,
243 NodeSet nodeset,
244 double* coeff_out,
245 size_t* indices_out,
246 size_t& num_coeff,
247 MsqError& err ) const
248 {
249 if (nodeset.have_any_mid_node()) {
250 MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
251 return;
252 }
253
254 switch (loc.dimension) {
255 case 0:
256 coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
257 break;
258 case 1:
259 coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff );
260 break;
261 case 2:
262 coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff );
263 break;
264 case 3:
265 coefficients_at_mid_elem( coeff_out, indices_out, num_coeff );
266 break;
267 default:
268 MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
269 }
270 }
271
derivatives(Sample loc,NodeSet nodeset,size_t * vertex_indices_out,MsqVector<3> * d_coeff_d_xi_out,size_t & num_vtx,MsqError & err) const272 void LinearPyramid::derivatives( Sample loc,
273 NodeSet nodeset,
274 size_t* vertex_indices_out,
275 MsqVector<3>* d_coeff_d_xi_out,
276 size_t& num_vtx,
277 MsqError& err ) const
278 {
279 if (nodeset.have_any_mid_node()) {
280 MSQ_SETERR(err)(nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
281 return;
282 }
283
284 switch (loc.dimension) {
285 case 0:
286 if (loc.number == 4) {
287 set_equal_derivatives( 0.0, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
288 }
289 else {
290 set_corner_derivatives( loc.number, 1.0, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
291 }
292 break;
293 case 1:
294 if (loc.number < 4) {
295 set_edge_derivatives( loc.number, 0.50, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
296 }
297 else {
298 set_corner_derivatives( loc.number-4, 0.50, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
299 }
300 break;
301 case 2:
302 if (loc.number == 4) {
303 set_equal_derivatives( 0.5, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
304 }
305 else {
306 set_edge_derivatives( loc.number, 0.25, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
307 }
308 break;
309 case 3:
310 set_equal_derivatives( 0.25, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
311 break;
312 default:
313 MSQ_SETERR(err)("Invalid/unsupported logical dimension",MsqError::INVALID_ARG);
314 }
315 }
316
ideal(Sample location,MsqMatrix<3,3> & J,MsqError &) const317 void LinearPyramid::ideal( Sample location,
318 MsqMatrix<3,3>& J,
319 MsqError& ) const
320 {
321 // For an ideal element with unit edge length at the base and unit
322 // height, the Jacobian matrix is:
323 // | 1-zeta 0 1/2 - xi |
324 // | 0 1-zeta 1/2 - eta |
325 // | 0 0 1 |
326 //
327 // The coefficient to produce a unit determinant
328 // is therefore (1-zeta)^(-2/3).
329 //
330 // Thus the unit-determinate ideal element Jacobian
331 // is, given alpha = (1-zeta)^(-1/3):
332 //
333 // | 1/alpha 0 alpha^2 (1/2 - xi) |
334 // | 0 1/alpha alpha^2 (1/2 - eta)|
335 // | 0 0 alpha^2 |
336 //
337 // There are only three zeta values of interest:
338 // zeta = 1 : the degenerate case
339 // zeta = 0 : both 1/alpha and alpha^2 are 1.0
340 // zeta = 1/2 : 1/alpha = 1/cbrt(2.0) and alpha^2 = 2*(1/alpha)
341
342 // special case for apex
343 if (location.dimension == 0 && location.number == 4) {
344 J = MsqMatrix<3,3>(0.0);
345 return;
346 }
347
348 // These are always zero
349 J(0,1) = J(1,0) = J(2,0) = J(2,1) = 0.0;
350
351 // Set diagonal terms and magnitude of terms in 3rd column based on zeta
352
353 // All of the zeta=0 locations
354 double f;
355 if ( location.dimension == 0 ||
356 (location.dimension == 1 && location.number < 4) ||
357 (location.dimension == 2 && location.number == 4)) {
358 J(0,0) = J(1,1) = J(2,2) = 1.0;
359 f = 0.5;
360 }
361 // all of the zeta=1/2 locations
362 else {
363 f = J(0,0) = J(1,1) = 0.79370052598409979;
364 J(2,2) = 2.0*f;
365 }
366
367 // Set terms in 3rd column based on xi,eta
368
369 // The xi = eta = 0.5 locations (mid-element in xi and eta)
370 if ( location.dimension == 3 ||
371 (location.dimension == 2 && location.number == 4)) {
372 J(0,2) = J(1,2) = 0.0;
373 }
374 // The corner locations
375 else if ( location.dimension == 0 ||
376 (location.dimension == 1 && location.number >= 4)) {
377 switch (location.number % 4) {
378 case 0: J(0,2) = f; J(1,2) = f; break;
379 case 1: J(0,2) = -f; J(1,2) = f; break;
380 case 2: J(0,2) = -f; J(1,2) = -f; break;
381 case 3: J(0,2) = f; J(1,2) = -f; break;
382 }
383 }
384 // The mid-edge locations
385 else {
386 switch (location.number) {
387 case 0: J(0,2) = 0; J(1,2) = f; break;
388 case 1: J(0,2) = -f; J(1,2) = 0; break;
389 case 2: J(0,2) = 0; J(1,2) = -f; break;
390 case 3: J(0,2) = f; J(1,2) = 0; break;
391 }
392 }
393 }
394
395 } // namespace MBMesquite
396