b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_klt_generic_shape.H
1//------------------------------------------------------------------------
2// b2element_klt_generic_shape.H --
3//
4// Base class for element "shapes".
5//
6// Copyright (c) 2017
7// SMR Engineering & Development SA
8// 2502 Bienne, Switzerland
9//
10// All Rights Reserved. Proprietary source code. The contents of
11// this file may not be disclosed to third parties, copied or
12// duplicated in any form, in whole or in part, without the prior
13// written permission of SMR.
14//------------------------------------------------------------------------
15
16#ifndef B2ELEMENT_KLT_GENERIC_SHAPE_H_
17#define B2ELEMENT_KLT_GENERIC_SHAPE_H_
18
19#include <algorithm>
20#include <cassert>
21
22#include "elements/klt_shells/b2element_klt_compute_a1a2.H"
23#include "elements/klt_shells/b2element_klt_compute_a3.H"
24#include "elements/klt_shells/b2element_klt_util.H"
25#include "elements/smr/b2element_continuum_integrate.H"
26#include "elements/smr/b2element_continuum_shape.H"
27#include "elements/smr/b2element_continuum_util.H"
28#include "model/b2domain.H"
29#include "model/b2element.H"
30#include "model_imp/b2hnode.H"
31#include "utils/b2tensor_calculus.H"
32
33namespace b2000::klt {
34
35const bool disp_interpolate_poly = false;
36
37class GenericShape {
38public:
39 GenericShape(const size_t num_nodes_, const int num_edges_)
40 : num_nodes(num_nodes_), num_edges(num_edges_), has_weight(false), has_weight_dof(false) {
41 nodes.resize(num_nodes);
42 node_dof_offset.resize(num_nodes);
43 node_dof_num.resize(num_nodes);
44 node_coor.resize(3, num_nodes);
45 node_weights.resize(num_nodes);
46 node_indices_all.reserve(num_nodes);
47 for (size_t i = 0; i != num_nodes; ++i) { node_indices_all.push_back(i); }
48 }
49
50 size_t get_num_nodes() { return num_nodes; }
51
52 int get_num_edges() { return num_edges; }
53
54 virtual IntegrationScheme create_ischeme() const = 0;
55
56 virtual int get_edge_order(const int edge) const = 0;
57
58 virtual int get_edge_shape_order(const int edge) const = 0;
59
60 virtual void get_xi_all_from_xi_edge(
61 const int edge, const double xi, double xi_all[2]) const = 0;
62
63 virtual IntegrationScheme create_ischeme_edge(const int edge) const = 0;
64
65 virtual void get_edge_direction(const int edge, double edge_dir[2]) const = 0;
66
67 virtual void get_node_indices(
68 const int sub_type, const int sub_index, b2linalg::Index& indices) const = 0;
69
70 virtual void compute_nodes_interpolation_straight_d0(
71 const double xi[2], b2linalg::Vector<double, b2linalg::Vdense>& N) const = 0;
72
73 virtual void compute_nodes_interpolation_straight_d1(
74 const double xi[2], b2linalg::Matrix<double, b2linalg::Mrectangle>& N) const = 0;
75
76 virtual void compute_nodes_interpolation_straight_d2(
77 const double xi[2], b2linalg::Matrix<double, b2linalg::Mrectangle>& N) const = 0;
78
79 virtual void compute_polynomial_basis_functions_internal(
80 const double xi[2], b2linalg::Matrix<double, b2linalg::Mrectangle>& N) const = 0;
81
82 void set_nodes(Element& element, std::pair<int, Node* const*> nodes_) {
83 if (size_t(nodes_.first) != num_nodes) {
84 SizeError() << "Element " << element.get_object_name() << " needs " << num_nodes
85 << " nodes." << THROW;
86 }
87 has_weight = false;
88 for (size_t i = 0; i != num_nodes; ++i) { nodes[i] = nodes_.second[i]; }
89 }
90
91 std::pair<int, Node* const*> get_nodes() const {
92 return std::pair<int, Node* const*>(num_nodes, &nodes[0]);
93 }
94
95 void compute_polynomial_basis_functions(
96 const b2linalg::Index& node_indices, const double xi[2], const int max_depth,
97 b2linalg::Matrix<double, b2linalg::Mrectangle>& basis_functions) const {
98 // Calculate all basis functions.
99 b2linalg::Matrix<double, b2linalg::Mrectangle> N_all;
100 compute_polynomial_basis_functions_internal(xi, N_all);
101
102 // Copy basis functions.
103 basis_functions.resize(node_indices.size(), max_depth);
104 for (int dd = 0; dd != max_depth; ++dd) {
105 for (size_t i = 0; i != node_indices.size(); ++i) {
106 basis_functions[dd][i] = N_all[dd][node_indices[i]];
107 }
108 }
109 }
110
111 void compute_polynomial_basis_functions_all(
112 const double xi[2],
113 b2linalg::Matrix<double, b2linalg::Mrectangle>& basis_functions) const {
114 basis_functions.resize(node_indices_all.size(), 6);
115 compute_polynomial_basis_functions_internal(xi, basis_functions);
116 }
117
118 virtual void compute_rational_basis_functions(
119 const b2linalg::Index& node_indices, const double xi[2], const int max_depth,
120 b2linalg::Matrix<double, b2linalg::Mrectangle>& basis_functions) {
121 b2linalg::Matrix<double, b2linalg::Mrectangle> P_all;
122 compute_polynomial_basis_functions_all(xi, P_all);
123
124 basis_functions.resize(node_indices.size(), max_depth);
125 if (!has_weight) {
126 for (size_t i = 0; i != node_indices.size(); ++i) {
127 const size_t ii = node_indices[i];
128 for (int dd = 0; dd != max_depth; ++dd) { basis_functions[dd][i] = P_all[dd][ii]; }
129 }
130 } else {
131 double W[6] = {};
132 for (int dd = 0; dd != 6; ++dd) {
133 W[dd] = 0;
134 for (size_t i = 0; i != num_nodes; ++i) { W[dd] += P_all[dd][i] * node_weights[i]; }
135 }
136 const double WI = 1. / W[d00];
137 const double WI2 = WI * WI;
138 const double WI3 = WI * WI * WI;
139
140 for (size_t i = 0; i != node_indices.size(); ++i) {
141 const size_t ii = node_indices[i];
142 double Q[6];
143 for (int dd = 0; dd != 6; ++dd) { Q[dd] = P_all[dd][ii] * node_weights[ii]; }
144 double R[6];
145
146 R[d00] = Q[d00] * WI;
147 R[d10] = Q[d10] * WI - Q[d00] * W[d10] * WI2;
148 R[d01] = Q[d01] * WI - Q[d00] * W[d01] * WI2;
149 R[d20] = Q[d20] * WI - Q[d10] * W[d10] * WI2 - Q[d10] * W[d10] * WI2
150 - Q[d00] * W[d20] * WI2 + 2 * Q[d00] * W[d10] * W[d10] * WI3;
151 R[d11] = Q[d11] * WI - Q[d10] * W[d01] * WI2 - Q[d01] * W[d10] * WI2
152 - Q[d00] * W[d11] * WI2 + 2 * Q[d00] * W[d10] * W[d01] * WI3;
153 R[d02] = Q[d02] * WI - Q[d01] * W[d01] * WI2 - Q[d01] * W[d01] * WI2
154 - Q[d00] * W[d02] * WI2 + 2 * Q[d00] * W[d01] * W[d01] * WI3;
155
156 for (int dd = 0; dd != max_depth; ++dd) { basis_functions[dd][i] = R[dd]; }
157 }
158 }
159 }
160
161 void init(Model& model) {
162 if (dof_type.empty()) {
163 has_weight = false;
164 for (size_t i = 0; i != num_nodes; ++i) {
165 // Get vertex coordinates in branch-global reference
166 // frame and get weights.
167 nodes[i]->get_coor(&node_coor[i][0]);
168 node_weights[i] = nodes[i]->get_weight();
169 has_weight = (has_weight || (node_weights[i] != 1));
170
171 // Get the number of dof and the dof types for each node.
172 int b[4] = {DOF_TYPE_DX, DOF_TYPE_DY, DOF_TYPE_DZ, DOF_TYPE_UNKNOWN};
173 node_dof_offset[i] = dof_numbering.size();
174 if (i >= size_t(num_edges) && NodeType4::is_dof_subset_s(nodes[i])) {
175 b[3] = DOF_TYPE_DW;
176 node_dof_num[i] = 4;
177 has_weight_dof = true;
178 } else {
179 assert(NodeType3::is_dof_subset_s(nodes[i]));
180 node_dof_num[i] = 3;
181 }
182 dof_type.insert(dof_type.end(), b, b + node_dof_num[i]);
183 }
184 }
185 }
186
187 void init_dof_numbering() {
188 assert(!dof_type.empty());
189
190 if (dof_numbering.empty()) {
191 dof_numbering.reserve(dof_type.size());
192 for (size_t i = 0; i != num_nodes; ++i) {
193 size_t a[4] = {};
194 node_dof_offset[i] = dof_numbering.size();
195 if (i >= size_t(num_edges) && NodeType4::is_dof_subset_s(nodes[i])) {
196 NodeType4::get_global_dof_numbering_s(a, nodes[i]);
197 } else {
198 assert(NodeType3::is_dof_subset_s(nodes[i]));
199 NodeType3::get_global_dof_numbering_s(a, nodes[i]);
200 }
201 dof_numbering.insert(dof_numbering.end(), a, a + node_dof_num[i]);
202 }
203 }
204 assert(dof_numbering.size() == dof_type.size());
205 }
206
207 void get_dof_info(
208 const b2linalg::Index& node_indices, b2linalg::Index& dof_numbering_,
209 std::vector<int>& dof_type_) {
210 init_dof_numbering();
211 dof_numbering_.clear();
212 dof_type_.clear();
213 for (size_t i = 0; i != node_indices.size(); ++i) {
214 const size_t ii = node_indices[i];
215 const size_t a = node_dof_offset[ii];
216 const size_t n = node_dof_num[ii];
217 dof_numbering_.insert(dof_numbering_.end(), &dof_numbering[a], &dof_numbering[a] + n);
218 dof_type_.insert(dof_type_.end(), &dof_type[a], &dof_type[a] + n);
219 }
220 }
221
222 virtual void interpolate_and_transform_midsurface_initial(
223 const b2linalg::Index& node_indices, const int geom_type, const double geom_param[2],
224 const int max_depth, const b2linalg::Matrix<double, b2linalg::Mrectangle>& N,
225 Vec3d Pos[6]) {
226 assert(max_depth >= 1 && max_depth <= 6);
227 assert(N.size1() == node_indices.size());
228 assert(N.size2() == size_t(max_depth));
229
230 // Interpolate initial and deformed midsurface and
231 // derivatives w.r.t. the parametric coordinates.
232 for (int dd = 0; dd != max_depth; ++dd) {
233 Pos[dd].set_zero();
234 for (size_t i = 0; i != node_indices.size(); ++i) {
235 const size_t ii = node_indices[i];
236 for (int j = 0; j != 3; ++j) { Pos[dd][j] += node_coor[ii][j] * N[dd][i]; }
237 }
238 }
239
240 //
241 // Transformation to shell.
242 //
243 if (geom_type == GEOM_ANALYTIC_CYLINDER) {
244 Vec3d F[6];
245 for (int dd = 0; dd != max_depth; ++dd) { F[dd] = Pos[dd]; }
246
247 const double radius = geom_param[0];
248 const double alpha_rad = F[d00][1] / radius;
249 const double sa = std::sin(alpha_rad);
250 const double ca = std::cos(alpha_rad);
251
252 const double rz[6] = {
253 radius - F[d00][2], -F[d10][2], -F[d01][2], -F[d20][2], -F[d11][2], -F[d02][2],
254 };
255
256 const double sad[6] = {
257 sa,
258 ca * F[d10][1] / radius,
259 ca * F[d01][1] / radius,
260 -sa * F[d10][1] * F[d10][1] / radius / radius + ca * F[d20][1] / radius,
261 -sa * F[d10][1] * F[d01][1] / radius / radius + ca * F[d11][1] / radius,
262 -sa * F[d01][1] * F[d01][1] / radius / radius + ca * F[d02][1] / radius,
263 };
264
265 const double cad[6] = {
266 ca,
267 -sa * F[d10][1] / radius,
268 -sa * F[d01][1] / radius,
269 -ca * F[d10][1] * F[d10][1] / radius / radius - sa * F[d20][1] / radius,
270 -ca * F[d10][1] * F[d01][1] / radius / radius - sa * F[d11][1] / radius,
271 -ca * F[d01][1] * F[d01][1] / radius / radius - sa * F[d02][1] / radius,
272 };
273
274 Pos[d00][0] = F[d00][0];
275 Pos[d00][1] = rz[d00] * sad[d00];
276 Pos[d00][2] = radius - rz[d00] * cad[d00];
277
278 if (max_depth >= d10) {
279 Pos[d10][0] = F[d10][0];
280 Pos[d10][1] = rz[d10] * sad[d00] + rz[d00] * sad[d10];
281 Pos[d10][2] = -rz[d10] * cad[d00] - rz[d00] * cad[d10];
282 }
283
284 if (max_depth >= d01) {
285 Pos[d01][0] = F[d01][0];
286 Pos[d01][1] = rz[d01] * sad[d00] + rz[d00] * sad[d01];
287 Pos[d01][2] = -rz[d01] * cad[d00] - rz[d00] * cad[d01];
288 }
289
290 if (max_depth >= d20) {
291 Pos[d20][0] = F[d20][0];
292 Pos[d20][1] = rz[d20] * sad[d00] + rz[d10] * sad[d10] + rz[d10] * sad[d10]
293 + rz[d00] * sad[d20];
294 Pos[d20][2] = -rz[d20] * cad[d00] - rz[d10] * cad[d10] - rz[d10] * cad[d10]
295 - rz[d00] * cad[d20];
296 }
297
298 if (max_depth >= d11) {
299 Pos[d11][0] = F[d11][0];
300 Pos[d11][1] = rz[d11] * sad[d00] + rz[d10] * sad[d01] + rz[d01] * sad[d10]
301 + rz[d00] * sad[d11];
302 Pos[d11][2] = -rz[d11] * cad[d00] - rz[d10] * cad[d01] - rz[d01] * cad[d10]
303 - rz[d00] * cad[d11];
304 }
305
306 if (max_depth >= d02) {
307 Pos[d02][0] = F[d02][0];
308 Pos[d02][1] = rz[d02] * sad[d00] + rz[d10] * sad[d10] + rz[d10] * sad[d10]
309 + rz[d00] * sad[d02];
310 Pos[d02][2] = -rz[d02] * cad[d00] - rz[d10] * cad[d10] - rz[d10] * cad[d10]
311 - rz[d00] * cad[d02];
312 }
313 } else if (geom_type == GEOM_ANALYTIC_HYPERBOLIC_PARABOLOID) {
314 Vec3d F[6];
315 for (int dd = 0; dd != max_depth; ++dd) { F[dd] = Pos[dd]; }
316
317 const double x = F[d00][0];
318 const double y = F[d00][1];
319 const double z = F[d00][2];
320 assert(z == 0);
321
322 const Vec3d& Fr = F[d10];
323 const Vec3d& Fs = F[d01];
324
325 Pos[d00][0] = x;
326 Pos[d00][1] = y;
327 Pos[d00][2] = x * x - y * y;
328
329 if (max_depth >= d10) {
330 Pos[d10][0] = Fr[0];
331 Pos[d10][1] = Fr[1];
332 Pos[d10][2] = 2 * x * Fr[0] - 2 * y * Fr[1];
333 }
334
335 if (max_depth >= d01) {
336 Pos[d01][0] = Fs[0];
337 Pos[d01][1] = Fs[1];
338 Pos[d01][2] = 2 * x * Fs[0] - 2 * y * Fs[1];
339 }
340
341 if (max_depth >= d20) {
342 Pos[d20][0] = 0;
343 Pos[d20][1] = 0;
344 Pos[d20][2] = 2 * Fr[0] * Fr[0] - 2 * Fr[1] * Fr[1];
345 }
346
347 if (max_depth >= d11) {
348 Pos[d11][0] = 0;
349 Pos[d11][1] = 0;
350 Pos[d11][2] = 2 * Fr[0] * Fs[0] - 2 * Fr[1] * Fs[1];
351 }
352
353 if (max_depth >= d02) {
354 Pos[d02][0] = 0;
355 Pos[d02][1] = 0;
356 Pos[d02][2] = 2 * Fs[0] * Fs[0] - 2 * Fs[1] * Fs[1];
357 }
358 }
359 }
360
361 void interpolate_and_transform_midsurface_initial_and_deformed(
362 const b2linalg::Index& node_indices, const int geom_type, const double geom_param[2],
363 const int max_depth, const b2linalg::Matrix<double, b2linalg::Mrectangle>& N,
364 Vec3d Pos[6], const b2linalg::Index& dof_numbering,
365 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const bool linear,
366 XiDofVec3d<6>& pos) {
367 interpolate_and_transform_midsurface_initial(
368 node_indices, geom_type, geom_param, max_depth, N, Pos);
369
370 assert(dof_numbering.size() == 3 * node_indices.size());
371
372 b2linalg::Index node_dof_offset_red;
373 b2linalg::Index node_dof_num_red;
374 {
375 size_t offset = 0;
376 for (size_t i = 0; i != node_indices.size(); ++i) {
377 const size_t ii = node_indices[i];
378 node_dof_offset_red.push_back(offset);
379 node_dof_num_red.push_back(node_dof_num[ii]);
380 offset += node_dof_num[ii];
381 }
382 assert(offset = dof_numbering.size());
383 }
384
385 pos.resize(dof_numbering.size(), true, false);
386 for (int dd = 0; dd != max_depth; ++dd) {
387 pos[dd].d0 = Pos[dd];
388 for (size_t i = 0; i != node_indices.size(); ++i) {
389 const size_t a = node_dof_offset_red[i];
390 const size_t b = a + node_dof_num_red[i];
391
392 assert(node_dof_num_red[i] == 3 && a % 3 == 0);
393 for (size_t j = a; j != b; ++j) {
394 pos[dd].d1[j][j % 3] = N[dd][i];
395 if (!linear && dof.size2() > 0) {
396 pos[dd].d0[j % 3] += N[dd][i] * dof[0][dof_numbering[j]];
397 }
398 }
399 }
400 }
401 }
402
403 virtual void compute_midsurface_initial_new_d0(
404 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
405 const double geom_param[2], Vec3d& Pos,
406 b2linalg::Matrix<double, b2linalg::Mrectangle>& N) {
407 compute_rational_basis_functions(node_indices, xi, 1, N);
408 interpolate_and_transform_midsurface_initial(
409 node_indices, geom_type, geom_param, 1, N, &Pos);
410
411 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 1, N); }
412 }
413
414 virtual void compute_midsurface_initial_new_d1(
415 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
416 const double geom_param[2], Vec3d Pos[3],
417 b2linalg::Matrix<double, b2linalg::Mrectangle>& N) {
418 compute_rational_basis_functions(node_indices, xi, 3, N);
419 interpolate_and_transform_midsurface_initial(
420 node_indices, geom_type, geom_param, 3, N, Pos);
421
422 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 3, N); }
423 }
424
425 virtual void compute_midsurface_initial_new_d2(
426 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
427 const double geom_param[2], Vec3d Pos[6],
428 b2linalg::Matrix<double, b2linalg::Mrectangle>& N) {
429 compute_rational_basis_functions(node_indices, xi, 6, N);
430 interpolate_and_transform_midsurface_initial(
431 node_indices, geom_type, geom_param, 6, N, Pos);
432
433 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 6, N); }
434 }
435
436 void compute_midsurface_initial_d0(
437 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
438 const double geom_param[2], Vec3d& Pos) {
439 b2linalg::Matrix<double, b2linalg::Mrectangle> N;
440 compute_rational_basis_functions(node_indices, xi, 1, N);
441 Vec3d PosV[6];
442 interpolate_and_transform_midsurface_initial(
443 node_indices, geom_type, geom_param, 1, N, PosV);
444 Pos = PosV[d00];
445
446 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 1, N); }
447 }
448
449 void compute_midsurface_initial_d1(
450 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
451 const double geom_param[2], Vec3d& Pos, Vec3d& A1, Vec3d& A2, Vec3d& A3) {
452 b2linalg::Matrix<double, b2linalg::Mrectangle> N;
453 compute_rational_basis_functions(node_indices, xi, 3, N);
454 Vec3d PosV[6];
455 interpolate_and_transform_midsurface_initial(
456 node_indices, geom_type, geom_param, 3, N, PosV);
457 Pos = PosV[d00];
458 A1 = PosV[d10];
459 A2 = PosV[d01];
460 A3 = normalised(outer_product(A1, A2));
461
462 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 3, N); }
463 }
464
465 void compute_midsurface_initial_d2(
466 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
467 const double geom_param[2], Vec3d& Pos, Vec3d A1[3], Vec3d A2[3], Vec3d A3[3]) {
468 b2linalg::Matrix<double, b2linalg::Mrectangle> N;
469 compute_rational_basis_functions(node_indices, xi, 6, N);
470 Vec3d PosV[6];
471 interpolate_and_transform_midsurface_initial(
472 node_indices, geom_type, geom_param, 6, N, PosV);
473 Pos = PosV[d00];
474 A1[d00] = PosV[d10];
475 A1[d10] = PosV[d20];
476 A1[d01] = PosV[d11];
477 A2[d00] = PosV[d01];
478 A2[d10] = PosV[d11];
479 A2[d01] = PosV[d02];
480 compute_A3(
481 (const double(*)[3])A1[0].data(), (const double(*)[3])A2[0].data(),
482 (double(*)[3])A3[0].data());
483
484 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 6, N); }
485 }
486
487 void compute_midsurface_initial_and_deformed_d0(
488 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
489 const double geom_param[2],
490 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const bool linear,
491 b2linalg::Matrix<double, b2linalg::Mrectangle>& N, Vec3d& Pos, XiDofVec3d<1>& pos) {
492 compute_rational_basis_functions(node_indices, xi, 1, N);
493 Vec3d PosV[6];
494 XiDofVec3d<6> posv;
495 interpolate_and_transform_midsurface_initial_and_deformed(
496 node_indices, geom_type, geom_param, 1, N, PosV, dof_numbering, dof, linear, posv);
497 Pos = PosV[d00];
498
499 pos.resize(dof_numbering.size(), true, false);
500
501 pos[d00] = pos[d00];
502
503 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 1, N); }
504 }
505
506 void compute_midsurface_initial_and_deformed_d1(
507 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
508 const double geom_param[2],
509 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const bool linear,
510 b2linalg::Matrix<double, b2linalg::Mrectangle>& N, Vec3d& Pos, Vec3d& A1, Vec3d& A2,
511 XiDofVec3d<1>& pos, XiDofVec3d<1>& a1, XiDofVec3d<1>& a2) {
512 compute_rational_basis_functions(node_indices, xi, 3, N);
513 Vec3d PosV[6];
514 XiDofVec3d<6> posv;
515 interpolate_and_transform_midsurface_initial_and_deformed(
516 node_indices, geom_type, geom_param, 3, N, PosV, dof_numbering, dof, linear, posv);
517 Pos = PosV[d00];
518 A1 = PosV[d10];
519 A2 = PosV[d01];
520
521 pos.resize(dof_numbering.size(), true, false);
522 a1.resize(dof_numbering.size(), true, false);
523 a2.resize(dof_numbering.size(), true, false);
524
525 pos[d00] = posv[d00];
526 a1[d00] = posv[d10];
527 a2[d00] = posv[d01];
528
529 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 3, N); }
530 }
531
532 void compute_midsurface_initial_and_deformed_d2(
533 const double xi[2], const b2linalg::Index& node_indices, const int geom_type,
534 const double geom_param[2],
535 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, const bool linear,
536 b2linalg::Matrix<double, b2linalg::Mrectangle>& N, Vec3d& Pos, Vec3d A1[3], Vec3d A2[3],
537 XiDofVec3d<1>& pos, XiDofVec3d<3>& a1, XiDofVec3d<3>& a2) {
538 compute_rational_basis_functions(node_indices, xi, 6, N);
539 Vec3d PosV[6];
540 XiDofVec3d<6> posv;
541 interpolate_and_transform_midsurface_initial_and_deformed(
542 node_indices, geom_type, geom_param, 6, N, PosV, dof_numbering, dof, linear, posv);
543 Pos = PosV[d00];
544 A1[d00] = PosV[d10];
545 A1[d10] = PosV[d20];
546 A1[d01] = PosV[d11];
547 A2[d00] = PosV[d01];
548 A2[d10] = PosV[d11];
549 A2[d01] = PosV[d02];
550
551 pos.resize(dof_numbering.size(), true, false);
552 a1.resize(dof_numbering.size(), true, false);
553 a2.resize(dof_numbering.size(), true, false);
554
555 pos[d00] = posv[d00];
556 a1[d00] = posv[d10];
557 a1[d10] = posv[d20];
558 a1[d01] = posv[d11];
559 a2[d00] = posv[d01];
560 a2[d10] = posv[d11];
561 a2[d01] = posv[d02];
562
563 if (disp_interpolate_poly) { compute_polynomial_basis_functions(node_indices, xi, 6, N); }
564 }
565
566public:
567 typedef node::HNode<
568 coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z> >,
569 coordof::Dof<coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ> > >
570 NodeType3;
571 typedef node::HNode<
572 coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z> >,
573 coordof::Dof<coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>, coordof::DW> >
574 NodeType4;
575
576 size_t num_nodes;
577 int num_edges;
578 std::vector<Node*> nodes;
579 bool has_weight;
580 bool has_weight_dof;
581 b2linalg::Index node_dof_offset;
582 b2linalg::Index node_dof_num;
583 b2linalg::Matrix<double, b2linalg::Mrectangle> node_coor;
584 b2linalg::Vector<double, b2linalg::Vdense> node_weights;
585 std::vector<int> dof_type;
586 b2linalg::Index dof_numbering;
587 b2linalg::Index node_indices_all;
588};
589
590} // namespace b2000::klt
591
592#endif // B2ELEMENT_KLT_GENERIC_SHAPE_H_
THROW
#define THROW
Definition b2exception.H:198
b2tensor_calculus.H
b2000::SizeError
GenericException< SizeError_name > SizeError
Definition b2exception.H:346