suanPan
shape.h
Go to the documentation of this file.
1/*******************************************************************************
2 * Copyright (C) 2017-2024 Theodore Chang
3 *
4 * This program is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program. If not, see <http://www.gnu.org/licenses/>.
16 ******************************************************************************/
29#ifndef SHAPE_H
30#define SHAPE_H
31
32#include <suanPan.h>
33
34namespace interpolation {
35 template<typename T> Row<T> linear(const T X, const T Y) { return {1., X, Y, X * Y}; }
36
37 template<typename T> Row<T> linear(const T X, const T Y, const T Z) { return {1., X, Y, Z, X * Y, Y * Z, Z * X}; }
38
39 template<typename T> Row<T> linear(const Col<T>& C) { return C.size() == 2 ? linear(C(0), C(1)) : linear(C(0), C(1), C(2)); }
40
41 template<typename T> Row<T> quadratic(const T X, const T Y) { return {1., X, Y, X * X, X * Y, Y * Y, X * X * Y, X * Y * Y, X * X * Y * Y}; }
42
43 template<typename T> Row<T> quadratic(const Col<T>& C) { return quadratic(C(0), C(1)); }
44
45 template<typename T> Row<T> cubic(const T X, const T Y) { return {1., X, Y, X * X, X * Y, Y * Y, X * X * X, X * X * Y, X * Y * Y, Y * Y * Y}; }
46
47 template<typename T> Row<T> cubic(const Col<T>& C) { return cubic(C(0), C(1)); }
48} // namespace interpolation
49
50namespace area {
51 template<typename T> T triangle(const Mat<T>& EC);
52
53 template<typename T> T shoelace(const Mat<T>& C);
54} // namespace area
55
56namespace shape {
57 template<typename T> Mat<T> truss(T int_pts, unsigned order = 0, unsigned num_node = 2);
58 template<typename T> Col<T> beam(T int_pts, unsigned order, double length);
59 template<typename T> Mat<T> triangle(const Col<T>& int_pts, unsigned order);
60 template<typename T> Mat<T> quad(const Mat<T>& int_pts, unsigned order, unsigned num_node = 4);
61 template<typename T> Mat<T> cube(const Mat<T>& int_pts, unsigned order, unsigned num_node = 8);
62
63 namespace plate {
64 template<typename T> Mat<T> triangle(const Col<T>& int_pts, unsigned order, unsigned num_node, const Mat<T>& nodes);
65 template<typename T> Mat<T> quad(const Col<T>& int_pts, unsigned order, unsigned num_node = 4);
66 } // namespace plate
67
68 template<typename T> Mat<T> stress(T X, T Y, unsigned S);
69
70 template<typename T> Mat<T> stress(const Col<T>& C, unsigned S);
71 template<typename T> Mat<T> stress5(const Col<T>& C);
72 template<typename T> Mat<T> stress7(const Col<T>& C);
73 template<typename T> Mat<T> stress9(const Col<T>& C);
74 template<typename T> Mat<T> stress11(const Col<T>& C);
75
76 template<typename T> Mat<T> stress5(T X, T Y);
77 template<typename T> Mat<T> stress7(T X, T Y);
78 template<typename T> Mat<T> stress9(T X, T Y);
79 template<typename T> Mat<T> stress11(T X, T Y);
80
81 inline mat stress5(const vec& C);
82 inline mat stress7(const vec& C);
83 inline mat stress9(const vec& C);
84 inline mat stress11(const vec& C);
85
86 template<typename T> Mat<T> strain(T X, T Y, T V, unsigned S);
87
88 template<typename T> Mat<T> strain(const Col<T>& C, T V, unsigned S);
89 template<typename T> Mat<T> strain5(T X, T Y, T V);
90 template<typename T> Mat<T> strain7(T X, T Y, T V);
91 template<typename T> Mat<T> strain9(T X, T Y, T V);
92 template<typename T> Mat<T> strain11(T X, T Y, T V);
93
94 template<typename T> Mat<T> strain5(const Col<T>& C, T V);
95 template<typename T> Mat<T> strain7(const Col<T>& C, T V);
96 template<typename T> Mat<T> strain9(const Col<T>& C, T V);
97 template<typename T> Mat<T> strain11(const Col<T>& C, T V);
98
99 inline mat strain5(const vec& C, double V);
100 inline mat strain7(const vec& C, double V);
101 inline mat strain9(const vec& C, double V);
102 inline mat strain11(const vec& C, double V);
103
104 template<typename T> Mat<T> linear_stress(T X, T Y);
105} // namespace shape
106
107template<typename T> T area::triangle(const Mat<T>& EC) { return .5 * (EC(0, 0) * (EC(1, 1) - EC(2, 1)) + EC(1, 0) * (EC(2, 1) - EC(0, 1)) + EC(2, 0) * (EC(0, 1) - EC(1, 1))); }
108
109template<typename T> T area::shoelace(const Mat<T>& C) {
110 suanpan_assert([&] { if(2 != C.n_cols) throw invalid_argument("need two columns"); });
111
112 const auto S = C.n_rows;
113 Mat<T> E = arma::resize(C, S + 1, C.n_cols);
114
115 E.tail_rows(1) = C.head_rows(1);
116
117 const vec& X = E.col(0);
118 const vec& Y = E.col(1);
119
120 return .5 * fabs(arma::dot(X.head(S), Y.tail(S)) - arma::dot(X.tail(S), Y.head(S)));
121}
122
123template<typename T> Mat<T> shape::truss(const T int_pts, const unsigned order, const unsigned num_node) {
124 Mat<T> N(1, num_node);
125
126 if(const auto& X = int_pts; num_node == 2) {
127 if(order == 0) {
128 N(0, 0) = 1. - X;
129 N(0, 1) = 1. + X;
130 }
131 else N(0, 0) = -(N(0, 1) = 1.);
132 N *= .5;
133 }
134 else if(num_node == 3) {
135 if(order == 0) {
136 const auto XX = X * X;
137 N(0, 0) = .5 * (XX - X);
138 N(0, 1) = 1. - XX;
139 N(0, 2) = .5 * (XX + X);
140 }
141 else {
142 N(0, 0) = X - .5;
143 N(0, 1) = -2. * X;
144 N(0, 2) = X + .5;
145 }
146 }
147
148 return N;
149}
150
151template<typename T> Col<T> shape::beam(const T int_pts, const unsigned order, const double length) {
152 Col<T> N(4);
153
154 const auto XP = 1. + int_pts;
155 const auto XM = 1. - int_pts;
156 const auto XPP = XP * XP;
157 const auto XMM = XM * XM;
158
159 if(order == 0) {
160 N(0) = 2. * XMM * (XP + 1.);
161 N(1) = length * XMM * XP;
162 N(2) = 2. * XPP * (XM + 1.);
163 N(3) = length * XM * XPP;
164 }
165 else if(order == 1) {
166 N(0) = -6. * XP * XM;
167 N(1) = length * XM * (3. * int_pts + 1.);
168 N(2) = 6. * XP * XM;
169 N(3) = length * XP * (3. * int_pts - 1.);
170 }
171
172 N *= .125;
173
174 return N;
175}
176
184template<typename T> Mat<T> shape::triangle(const Col<T>& int_pts, const unsigned order) {
185 Mat<T> N;
186
187 if(order != 0 && order != 1) throw invalid_argument("order needs to be either 0 or 1");
188
189 N.zeros(order + 1llu, 6);
190
191 if(const auto &X = int_pts(0), &Y = int_pts(1); order == 0) {
192 N(0, 0) = 1.;
193 N(0, 1) = X;
194 N(0, 2) = Y;
195 N(0, 3) = X * Y;
196 N(0, 4) = X * X;
197 N(0, 5) = Y * Y;
198 }
199 else if(order == 1) {
200 N(0, 1) = N(1, 2) = 1.;
201 N(0, 4) = 2. * (N(1, 3) = X);
202 N(1, 5) = 2. * (N(0, 3) = Y);
203 }
204
205 return N;
206}
207
208template<typename T> Mat<T> shape::quad(const Mat<T>& int_pts, const unsigned order, const unsigned num_node) {
209 Mat<T> N;
210
211 if(order != 0 && order != 1) throw invalid_argument("order needs to be either 0 or 1");
212 if(num_node < 4 || num_node > 8) throw invalid_argument("number of nodes must between 4 and 8");
213
214 N.zeros(order + 1llu, num_node);
215
216 const auto& X = int_pts(0);
217 const auto& Y = int_pts(1);
218
219 if(const auto XP = 1. + X, XM = 1. - X, YP = 1. + Y, YM = 1. - Y; 8 == num_node) {
220 if(const auto XX = X * X, YY = Y * Y, XY = X * Y; 0 == order) {
221 N(0, 7) = .5 * XM * (1. - YY);
222 N(0, 6) = .5 * (1. - XX) * YP;
223 N(0, 5) = .5 * XP * (1. - YY);
224 N(0, 4) = .5 * (1. - XX) * YM;
225 N(0, 0) = .25 * XM * YM - .5 * (N(0, 4) + N(0, 7));
226 N(0, 1) = .25 * XP * YM - .5 * (N(0, 4) + N(0, 5));
227 N(0, 2) = .25 * XP * YP - .5 * (N(0, 5) + N(0, 6));
228 N(0, 3) = .25 * XM * YP - .5 * (N(0, 6) + N(0, 7));
229 }
230 else if(1 == order) {
231 const auto X2 = .5 * X;
232 const auto Y2 = .5 * Y;
233 const auto X4 = .25 * X;
234 const auto Y4 = .25 * Y;
235 const auto X24 = .25 * XX;
236 const auto Y24 = .25 * YY;
237 const auto XY2 = .5 * XY;
238 N(1, 7) = XY - Y;
239 N(1, 6) = .5 - .5 * XX;
240 N(1, 5) = -Y - XY;
241 N(1, 4) = .5 * XX - .5;
242 N(1, 3) = Y2 - X4 - XY2 + X24;
243 N(1, 2) = X4 + Y2 + XY2 + X24;
244 N(1, 1) = Y2 - X4 + XY2 - X24;
245 N(1, 0) = X4 + Y2 - XY2 - X24;
246 N(0, 7) = .5 * YY - .5;
247 N(0, 6) = -X - XY;
248 N(0, 5) = .5 - .5 * YY;
249 N(0, 4) = XY - X;
250 N(0, 3) = X2 - Y4 + XY2 - Y24;
251 N(0, 2) = X2 + Y4 + XY2 + Y24;
252 N(0, 1) = X2 - Y4 - XY2 + Y24;
253 N(0, 0) = X2 + Y4 - XY2 - Y24;
254 }
255 }
256 else {
257 if(0 == order) {
258 N(0, 3) = XM * YP;
259 N(0, 2) = XP * YP;
260 N(0, 1) = XP * YM;
261 N(0, 0) = XM * YM;
262 }
263 else if(1 == order) {
264 N(1, 1) = -(N(1, 2) = XP);
265 N(1, 0) = -(N(1, 3) = XM);
266 N(0, 3) = -(N(0, 2) = YP);
267 N(0, 0) = -(N(0, 1) = YM);
268 }
269 N *= .25;
270 if(5 <= num_node) {
271 if(0 == order) {
272 N(0, 4) = .25 * (1. - X * X) * (1. - Y);
273 N(0, 0) -= .5 * N(0, 4);
274 N(0, 1) -= .5 * N(0, 4);
275 }
276 else {
277 N(0, 4) = -.5 * X * (1. - Y);
278 N(1, 4) = -.25 * (1. - X * X);
279 N(0, 0) -= .5 * N(0, 4);
280 N(0, 1) -= .5 * N(0, 4);
281 N(1, 0) -= .5 * N(1, 4);
282 N(1, 1) -= .5 * N(1, 4);
283 }
284 }
285 if(6 <= num_node) {
286 if(0 == order) {
287 N(0, 5) = .25 * (1. - Y * Y) * (1. + X);
288 N(0, 1) -= .5 * N(0, 5);
289 N(0, 2) -= .5 * N(0, 5);
290 }
291 else {
292 N(0, 5) = .25 * (1. - Y * Y);
293 N(1, 5) = -.5 * Y * (1. + X);
294 N(0, 1) -= .5 * N(0, 5);
295 N(0, 2) -= .5 * N(0, 5);
296 N(1, 1) -= .5 * N(1, 5);
297 N(1, 2) -= .5 * N(1, 5);
298 }
299 }
300 if(7 <= num_node) {
301 if(0 == order) {
302 N(0, 6) = .25 * (1. - X * X) * (1. + Y);
303 N(0, 2) -= .5 * N(0, 6);
304 N(0, 3) -= .5 * N(0, 6);
305 }
306 else {
307 N(0, 6) = -.5 * X * (1. + Y);
308 N(1, 6) = .25 * (1. - X * X);
309 N(0, 2) -= .5 * N(0, 6);
310 N(0, 3) -= .5 * N(0, 6);
311 N(1, 2) -= .5 * N(1, 6);
312 N(1, 3) -= .5 * N(1, 6);
313 }
314 }
315 }
316
317 return N;
318}
319
320template<typename T> Mat<T> shape::cube(const Mat<T>& int_pts, const unsigned order, const unsigned num_node) {
321 Mat<T> N;
322
323 if(order == 0) N.zeros(1, num_node);
324 else if(order == 1) N.zeros(3, num_node);
325 else throw invalid_argument("order needs to be either 0 or 1");
326
327 const auto& X = int_pts(0);
328 const auto& Y = int_pts(1);
329 const auto& Z = int_pts(2);
330
331 const auto XP = 1. + X;
332 const auto XM = 1. - X;
333 const auto YP = 1. + Y;
334 const auto YM = 1. - Y;
335 const auto ZP = 1. + Z;
336 const auto ZM = 1. - Z;
337
338 if(num_node == 8) {
339 if(order == 0) {
340 N(0, 0) = XM * YM * ZM;
341 N(0, 1) = XP * YM * ZM;
342 N(0, 2) = XP * YP * ZM;
343 N(0, 3) = XM * YP * ZM;
344 N(0, 4) = XM * YM * ZP;
345 N(0, 5) = XP * YM * ZP;
346 N(0, 6) = XP * YP * ZP;
347 N(0, 7) = XM * YP * ZP;
348 }
349 else if(order == 1) {
350 N(0, 0) = -YM * ZM;
351 N(0, 1) = YM * ZM;
352 N(0, 2) = YP * ZM;
353 N(0, 3) = -YP * ZM;
354 N(0, 4) = -YM * ZP;
355 N(0, 5) = YM * ZP;
356 N(0, 6) = YP * ZP;
357 N(0, 7) = -YP * ZP;
358 N(1, 0) = -XM * ZM;
359 N(1, 1) = -XP * ZM;
360 N(1, 2) = XP * ZM;
361 N(1, 3) = XM * ZM;
362 N(1, 4) = -XM * ZP;
363 N(1, 5) = -XP * ZP;
364 N(1, 6) = XP * ZP;
365 N(1, 7) = XM * ZP;
366 N(2, 0) = -XM * YM;
367 N(2, 1) = -XP * YM;
368 N(2, 2) = -XP * YP;
369 N(2, 3) = -XM * YP;
370 N(2, 4) = XM * YM;
371 N(2, 5) = XP * YM;
372 N(2, 6) = XP * YP;
373 N(2, 7) = XM * YP;
374 }
375 N *= .125;
376
377 return N;
378 }
379
380 if(num_node == 20) {
381 if(const auto XX = XP * XM, YY = YP * YM, ZZ = ZP * ZM; order == 0) {
382 N(0, 0) = .125 * XM * YM * ZM * (-2. - X - Y - Z);
383 N(0, 1) = .125 * XP * YM * ZM * (-2. + X - Y - Z);
384 N(0, 2) = .125 * XP * YP * ZM * (-2. + X + Y - Z);
385 N(0, 3) = .125 * XM * YP * ZM * (-2. - X + Y - Z);
386 N(0, 4) = .125 * XM * YM * ZP * (-2. - X - Y + Z);
387 N(0, 5) = .125 * XP * YM * ZP * (-2. + X - Y + Z);
388 N(0, 6) = .125 * XP * YP * ZP * (-2. + X + Y + Z);
389 N(0, 7) = .125 * XM * YP * ZP * (-2. - X + Y + Z);
390 N(0, 8) = .25 * XX * YM * ZM;
391 N(0, 9) = .25 * YY * XP * ZM;
392 N(0, 10) = .25 * XX * YP * ZM;
393 N(0, 11) = .25 * YY * XM * ZM;
394 N(0, 12) = .25 * XX * YM * ZP;
395 N(0, 13) = .25 * YY * XP * ZP;
396 N(0, 14) = .25 * XX * YP * ZP;
397 N(0, 15) = .25 * YY * XM * ZP;
398 N(0, 16) = .25 * ZZ * XM * YM;
399 N(0, 17) = .25 * ZZ * XP * YM;
400 N(0, 18) = .25 * ZZ * XP * YP;
401 N(0, 19) = .25 * ZZ * XM * YP;
402 }
403 else if(order == 1) {
404 N(0, 0) = YM * ZM * (2. * X + Y + Z + 1.) * .125;
405 N(0, 1) = YM * ZM * (2. * X - Y - Z - 1.) * .125;
406 N(0, 2) = YP * ZM * (2. * X + Y - Z - 1.) * .125;
407 N(0, 3) = YP * ZM * (2. * X - Y + Z + 1.) * .125;
408 N(0, 4) = YM * ZP * (2. * X + Y - Z + 1.) * .125;
409 N(0, 5) = YM * ZP * (2. * X - Y + Z - 1.) * .125;
410 N(0, 6) = YP * ZP * (2. * X + Y + Z - 1.) * .125;
411 N(0, 7) = YP * ZP * (2. * X - Y - Z + 1.) * .125;
412 N(1, 0) = XM * ZM * (2. * Y + X + Z + 1.) * .125;
413 N(1, 1) = XP * ZM * (2. * Y - X + Z + 1.) * .125;
414 N(1, 2) = XP * ZM * (2. * Y + X - Z - 1.) * .125;
415 N(1, 3) = XM * ZM * (2. * Y - X - Z - 1.) * .125;
416 N(1, 4) = XM * ZP * (2. * Y + X - Z + 1.) * .125;
417 N(1, 5) = XP * ZP * (2. * Y - X - Z + 1.) * .125;
418 N(1, 6) = XP * ZP * (2. * Y + X + Z - 1.) * .125;
419 N(1, 7) = XM * ZP * (2. * Y - X + Z - 1.) * .125;
420 N(2, 0) = XM * YM * (2. * Z + X + Y + 1.) * .125;
421 N(2, 1) = XP * YM * (2. * Z + Y - X + 1.) * .125;
422 N(2, 2) = XP * YP * (2. * Z - X - Y + 1.) * .125;
423 N(2, 3) = XM * YP * (2. * Z + X - Y + 1.) * .125;
424 N(2, 4) = XM * YM * (2. * Z - X - Y - 1.) * .125;
425 N(2, 5) = XP * YM * (2. * Z + X - Y - 1.) * .125;
426 N(2, 6) = XP * YP * (2. * Z + X + Y - 1.) * .125;
427 N(2, 7) = XM * YP * (2. * Z - X + Y - 1.) * .125;
428
429 N(0, 8) = -X * YM * ZM * .5;
430 N(0, 9) = YY * ZM * .25;
431 N(0, 10) = -X * YP * ZM * .5;
432 N(0, 11) = -YY * ZM * .25;
433 N(1, 8) = -XX * ZM * .25;
434 N(1, 9) = -Y * XP * ZM * .5;
435 N(1, 10) = XX * ZM * .25;
436 N(1, 11) = -Y * XM * ZM * .5;
437 N(2, 8) = -XX * YM * .25;
438 N(2, 9) = -XP * YY * .25;
439 N(2, 10) = -XX * YP * .25;
440 N(2, 11) = -XM * YY * .25;
441
442 N(0, 12) = -X * YM * ZP * .5;
443 N(0, 13) = YY * ZP * .25;
444 N(0, 14) = -X * YP * ZP * .5;
445 N(0, 15) = -YY * ZP * .25;
446 N(1, 12) = -XX * ZP * .25;
447 N(1, 13) = -Y * XP * ZP * .5;
448 N(1, 14) = XX * ZP * .25;
449 N(1, 15) = -Y * XM * ZP * .5;
450 N(2, 12) = XX * YM * .25;
451 N(2, 13) = XP * YY * .25;
452 N(2, 14) = XX * YP * .25;
453 N(2, 15) = XM * YY * .25;
454
455 N(0, 16) = -YM * ZZ * .25;
456 N(0, 17) = YM * ZZ * .25;
457 N(0, 18) = YP * ZZ * .25;
458 N(0, 19) = -YP * ZZ * .25;
459 N(1, 16) = -XM * ZZ * .25;
460 N(1, 17) = -XP * ZZ * .25;
461 N(1, 18) = XP * ZZ * .25;
462 N(1, 19) = XM * ZZ * .25;
463 N(2, 16) = -Z * XM * YM * .5;
464 N(2, 17) = -Z * XP * YM * .5;
465 N(2, 18) = -Z * XP * YP * .5;
466 N(2, 19) = -Z * XM * YP * .5;
467 }
468
469 return N;
470 }
471
472 throw invalid_argument("not supported");
473}
474
475template<typename T> Mat<T> shape::stress(const T X, const T Y, const unsigned S) {
476 Mat<T> N = zeros(3, S);
477
478 for(auto I = 0; I < 3; ++I) N(I, I) = 1.;
479
480 if(S >= 5) {
481 N(0, 4) = Y;
482 N(1, 3) = X;
483 if(S >= 7) {
484 N(2, 5) = -(N(0, 6) = X);
485 N(2, 6) = -(N(1, 5) = Y);
486 if(S >= 9) {
487 const auto X2 = X * X;
488 const auto Y2 = Y * Y;
489 const auto XY = X * Y;
490 N(1, 7) = N(0, 8) = 2. * XY;
491 N(2, 7) = -X2;
492 N(2, 8) = -Y2;
493 if(S == 11) {
494 N(1, 9) = 2. * X2 + (N(1, 10) = -Y2);
495 N(0, 10) = 2. * Y2 + (N(0, 9) = -X2);
496 N(2, 10) = N(2, 9) = 2. * XY;
497 }
498 }
499 }
500 }
501
502 return N;
503}
504
505template<typename T> Mat<T> shape::stress(const Col<T>& C, const unsigned S) { return stress(C(0), C(1), S); }
506
507template<typename T> Mat<T> shape::stress5(const Col<T>& C) { return stress(C, 5); }
508
509template<typename T> Mat<T> shape::stress7(const Col<T>& C) { return stress(C, 7); }
510
511template<typename T> Mat<T> shape::stress9(const Col<T>& C) { return stress(C, 9); }
512
513template<typename T> Mat<T> shape::stress11(const Col<T>& C) { return stress(C, 11); }
514
515template<typename T> Mat<T> shape::stress5(const T X, const T Y) { return stress(X, Y, 5); }
516
517template<typename T> Mat<T> shape::stress7(const T X, const T Y) { return stress(X, Y, 7); }
518
519template<typename T> Mat<T> shape::stress9(const T X, const T Y) { return stress(X, Y, 9); }
520
521template<typename T> Mat<T> shape::stress11(const T X, const T Y) { return stress(X, Y, 11); }
522
523template<typename T> Mat<T> shape::strain(const T X, const T Y, const T V, const unsigned S) {
524 Mat<T> N(3, S, fill::zeros);
525
526 N(0, 0) = N(1, 1) = 1.;
527
528 N(2, 2) = 2. + 2. * V;
529
530 N(0, 1) = N(1, 0) = -V;
531
532 if(S >= 5) {
533 N(0, 3) = -V * (N(1, 3) = X);
534 N(1, 4) = -V * (N(0, 4) = Y);
535 if(S >= 7) {
536 N(0, 5) = N(1, 4);
537 N(0, 6) = N(1, 3);
538
539 N(1, 5) = N(0, 4);
540 N(1, 6) = N(0, 3);
541
542 N(2, 5) = -X * N(2, 2);
543 N(2, 6) = -Y * N(2, 2);
544 if(S >= 9) {
545 const auto X2 = X * X, Y2 = Y * Y, XY = X * Y;
546
547 N(1, 8) = N(0, 7) = -V * (N(1, 7) = N(0, 8) = 2. * XY);
548
549 N(2, 7) = -X2 * N(2, 2);
550 N(2, 8) = -Y2 * N(2, 2);
551 if(S == 11) {
552 N(0, 9) = V * Y2 - (2. * V + 1.) * X2;
553 N(1, 9) = (2. + V) * X2 - Y2;
554
555 N(0, 10) = (2. + V) * Y2 - X2;
556 N(1, 10) = V * X2 - (2. * V + 1.) * Y2;
557
558 N(2, 10) = N(2, 9) = 2. * XY * N(2, 2);
559 }
560 }
561 }
562 }
563
564 return N;
565}
566
567template<typename T> Mat<T> shape::strain(const Col<T>& C, const T V, const unsigned S) { return strain(C(0), C(1), V, S); }
568
569template<typename T> Mat<T> shape::strain5(const T X, const T Y, const T V) { return strain(X, Y, V, 5); }
570
571template<typename T> Mat<T> shape::strain7(const T X, const T Y, const T V) { return strain(X, Y, V, 7); }
572
573template<typename T> Mat<T> shape::strain9(const T X, const T Y, const T V) { return strain(X, Y, V, 9); }
574
575template<typename T> Mat<T> shape::strain11(const T X, const T Y, const T V) { return strain(X, Y, V, 11); }
576
577template<typename T> Mat<T> shape::strain5(const Col<T>& C, const T V) { return strain(C, V, 5); }
578
579template<typename T> Mat<T> shape::strain7(const Col<T>& C, const T V) { return strain(C, V, 7); }
580
581template<typename T> Mat<T> shape::strain9(const Col<T>& C, const T V) { return strain(C, V, 9); }
582
583template<typename T> Mat<T> shape::strain11(const Col<T>& C, const T V) { return strain(C, V, 11); }
584
585template<typename T> Mat<T> shape::linear_stress(T X, T Y) {
586 Mat<T> N = eye(3, 9);
587 N.cols(3, 5) = X * eye(3, 3);
588 N.cols(6, 8) = Y * eye(3, 3);
589
590 return N;
591}
592
593template<typename T> Mat<T> shape::plate::triangle(const Col<T>& int_pts, const unsigned order, const unsigned num_node, const Mat<T>& nodes) {
594 suanpan_assert([&] { if(order > 1 || (num_node != 3 && num_node != 6) || (nodes.n_cols != 2) || nodes.n_rows < 3) throw invalid_argument("not supported"); });
595
596 Mat<T> N;
597
598 if(order == 0) N.zeros(1, 3llu * num_node);
599 else if(order == 1) N.zeros(3, 3llu * num_node);
600 else throw invalid_argument("order needs to be either 0 or 1");
601
602 Mat<T> TEMP(3, 3);
603 TEMP.row(0).fill(1.);
604 TEMP.rows(1, 2) = nodes.rows(0, 2).t();
605
606 const Mat<T> A = inv(TEMP);
607 const Col<T> W = solve(TEMP, vec{1., int_pts(0), int_pts(1)});
608
609 const vec L = std::initializer_list<double>{W(0), W(1), W(2), W(0), W(1)};
610 const vec B = std::initializer_list<double>{A(0, 1), A(1, 1), A(2, 1), A(0, 1), A(1, 1)};
611 const vec C = std::initializer_list<double>{A(0, 2), A(1, 2), A(2, 2), A(0, 2), A(1, 2)};
612
613 const auto DA = A.cols(1, 2);
614
615 if(3 == num_node) {
616 if(0 == order) {
617 auto IDX = 0;
618 for(auto I = 0; I < 3; ++I) {
619 const auto J = I + 1;
620 const auto K = J + 1;
621 N(IDX++) = L(I) * (1. - L(J) * L(J) - L(K) * L(K)) + L(I) * L(I) * (L(J) + L(K));
622 N(IDX++) = B(J) * L(I) * L(K) * (L(I) + .5 * L(J)) - B(K) * L(I) * L(J) * (L(I) + .5 * L(K));
623 N(IDX++) = C(J) * L(I) * L(K) * (L(I) + .5 * L(J)) - C(K) * L(I) * L(J) * (L(I) + .5 * L(K));
624 }
625 }
626 else {
627 mat DNDL(3, 3);
628 auto IDX = 0;
629 for(auto I = 0; I < 3; ++I) {
630 const auto J = I + 1;
631 const auto K = J + 1;
632
633 DNDL(0, 0) = L(J) + L(K);
634 DNDL(1, 1) = DNDL(2, 2) = -L(I);
635 DNDL(0, 1) = DNDL(1, 0) = L(I) - L(J);
636 DNDL(0, 2) = DNDL(2, 0) = L(I) - L(K);
637 DNDL(1, 2) = DNDL(2, 1) = 0.;
638
639 mat DD = DA.t() * DNDL * DA;
640
641 N(0, IDX) = -2. * DD(0, 0);
642 N(1, IDX) = -2. * DD(1, 1);
643 N(2, IDX++) = -2. * (DD(0, 1) + DD(1, 0));
644
645 DNDL(0, 0) = 4. * (B(J) * L(K) - B(K) * L(J));
646 DNDL(1, 1) = DNDL(2, 2) = 0.;
647 DNDL(0, 1) = DNDL(1, 0) = (B(J) - B(K)) * L(K) - 4 * B(K) * L(I);
648 DNDL(0, 2) = DNDL(2, 0) = 4. * B(J) * L(I) + (B(J) - B(K)) * L(J);
649 DNDL(1, 2) = DNDL(2, 1) = L(I) * (B(J) - B(K));
650
651 DD = DA.t() * DNDL * DA;
652
653 N(0, IDX) = -.5 * DD(0, 0);
654 N(1, IDX) = -.5 * DD(1, 1);
655 N(2, IDX++) = -.5 * (DD(0, 1) + DD(1, 0));
656
657 DNDL(0, 0) = 4. * (C(J) * L(K) - C(K) * L(J));
658 DNDL(0, 1) = DNDL(1, 0) = (C(J) - C(K)) * L(K) - 4 * C(K) * L(I);
659 DNDL(0, 2) = DNDL(2, 0) = 4. * C(J) * L(I) + (C(J) - C(K)) * L(J);
660 DNDL(1, 2) = DNDL(2, 1) = L(I) * (C(J) - C(K));
661
662 DD = DA.t() * DNDL * DA;
663
664 N(0, IDX) = -.5 * DD(0, 0);
665 N(1, IDX) = -.5 * DD(1, 1);
666 N(2, IDX++) = -.5 * (DD(0, 1) + DD(1, 0));
667 }
668 }
669 }
670
671 return N;
672}
673
674template<typename T> Mat<T> shape::plate::quad(const Col<T>& int_pts, const unsigned order, const unsigned num_node) {
675 Mat<T> N;
676
677 if(order == 0) N.zeros(1, 3llu * num_node);
678 else if(order == 1) N.zeros(2, 6llu * num_node);
679 else throw invalid_argument("order needs to be either 0 or 1");
680
681 const auto& X = int_pts(0);
682 const auto& Y = int_pts(1);
683
684 return N;
685}
686
687mat shape::stress5(const vec& C) { return stress(C, 5); }
688
689mat shape::stress7(const vec& C) { return stress(C, 7); }
690
691mat shape::stress9(const vec& C) { return stress(C, 9); }
692
693mat shape::stress11(const vec& C) { return stress(C, 11); }
694
695mat shape::strain5(const vec& C, const double V) { return strain(C, V, 5); }
696
697mat shape::strain7(const vec& C, const double V) { return strain(C, V, 7); }
698
699mat shape::strain9(const vec& C, const double V) { return strain(C, V, 9); }
700
701mat shape::strain11(const vec& C, const double V) { return strain(C, V, 11); }
702
703#endif
704
PlaneType::E
@ E
PlaneType::S
@ S
PlaneType::A
@ A
PlaneType::N
@ N
DOF::T
@ T
shape::strain9
Mat< T > strain9(T X, T Y, T V)
Definition: shape.h:573
shape::quad
Mat< T > quad(const Mat< T > &int_pts, unsigned order, unsigned num_node=4)
Definition: shape.h:208
shape::stress
Mat< T > stress(T X, T Y, unsigned S)
Definition: shape.h:475
shape::strain
Mat< T > strain(T X, T Y, T V, unsigned S)
Definition: shape.h:523
shape::strain
Mat< T > strain(const Col< T > &C, T V, unsigned S)
Definition: shape.h:567
shape::stress7
Mat< T > stress7(const Col< T > &C)
Definition: shape.h:509
shape::triangle
Mat< T > triangle(const Col< T > &int_pts, unsigned order)
compute the shape function or its derivative of six node triangle in global coordinate system
Definition: shape.h:184
shape::beam
Col< T > beam(T int_pts, unsigned order, double length)
Definition: shape.h:151
area::shoelace
T shoelace(const Mat< T > &C)
Definition: shape.h:109
area::triangle
T triangle(const Mat< T > &EC)
Definition: shape.h:107
shape::strain5
Mat< T > strain5(T X, T Y, T V)
Definition: shape.h:569
shape::truss
Mat< T > truss(T int_pts, unsigned order=0, unsigned num_node=2)
Definition: shape.h:123
shape::stress5
Mat< T > stress5(const Col< T > &C)
Definition: shape.h:507
shape::linear_stress
Mat< T > linear_stress(T X, T Y)
Definition: shape.h:585
shape::plate::quad
Mat< T > quad(const Col< T > &int_pts, unsigned order, unsigned num_node=4)
Definition: shape.h:674
shape::stress9
Mat< T > stress9(const Col< T > &C)
Definition: shape.h:511
shape::strain11
Mat< T > strain11(T X, T Y, T V)
Definition: shape.h:575
shape::cube
Mat< T > cube(const Mat< T > &int_pts, unsigned order, unsigned num_node=8)
Definition: shape.h:320
shape::stress11
Mat< T > stress11(const Col< T > &C)
Definition: shape.h:513
shape::stress
Mat< T > stress(const Col< T > &C, unsigned S)
Definition: shape.h:505
shape::plate::triangle
Mat< T > triangle(const Col< T > &int_pts, unsigned order, unsigned num_node, const Mat< T > &nodes)
Definition: shape.h:593
shape::strain7
Mat< T > strain7(T X, T Y, T V)
Definition: shape.h:571
area
Definition: shape.h:50
interpolation
Definition: shape.h:34
interpolation::quadratic
Row< T > quadratic(const T X, const T Y)
Definition: shape.h:41
interpolation::linear
Row< T > linear(const T X, const T Y)
Definition: shape.h:35
interpolation::cubic
Row< T > cubic(const T X, const T Y)
Definition: shape.h:45
shape
Definition: shape.h:56
suanpan_assert
void suanpan_assert(const std::function< void()> &F)
Definition: suanPan.h:296