b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_klt_compute_shape_q.H
1//------------------------------------------------------------------------
2// b2element_klt_compute_shape_q.H --
3//
4// Compute shape functions and derivatives for Bezier quadrilaterals.
5//
6// written by Thomas Ludwig
7// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
8//
9// (c) 2015,2016,2017 SMR Engineering & Development SA
10// 2502 Bienne, Switzerland
11//
12// (c) 2023 Deutsches Zentrum für Luft- und Raumfahrt (DLR) e.V.
13// Linder Höhe, 51147 Köln
14//
15// All Rights Reserved. Proprietary source code. The contents of
16// this file may not be disclosed to third parties, copied or
17// duplicated in any form, in whole or in part, without the prior
18// written permission of SMR.
19//------------------------------------------------------------------------
20
21#ifndef B2ELEMENT_KLT_COMPUTE_SHAPE_Q_H_
22#define B2ELEMENT_KLT_COMPUTE_SHAPE_Q_H_
23
24#include <algorithm>
25#include <cassert>
26
27#include "elements/klt_shells/b2element_klt_generic_shape.H"
28#include "elements/klt_shells/b2element_klt_util.H"
29#include "elements/smr/b2element_continuum_integrate.H"
30#include "elements/smr/b2element_continuum_shape.H"
31#include "elements/smr/b2element_continuum_util.H"
32#include "model/b2domain.H"
33#include "model/b2element.H"
34#include "model_imp/b2hnode.H"
35#include "utils/b2tensor_calculus.H"
36
37namespace b2000::klt {
38
39template <int ORDER_R, int ORDER_S, typename DATA>
40class QuadBezier : public GenericShape {
41public:
42 static const int num_nodes = (ORDER_R + 1) * (ORDER_S + 1);
43
44 static const int order_r = ORDER_R;
45
46 static const int order_s = ORDER_S;
47
48 static const size_t num_nodes_edge_r = ORDER_R + 1;
49
50 static const size_t num_nodes_edge_c1_r = (ORDER_R + 1) * 2;
51
52 static const size_t num_nodes_edge_c2_r = (ORDER_R + 1) * 3;
53
54 static const size_t num_nodes_edge_s = ORDER_S + 1;
55
56 static const size_t num_nodes_edge_c1_s = (ORDER_S + 1) * 2;
57
58 static const size_t num_nodes_edge_c2_s = (ORDER_S + 1) * 3;
59
60 QuadBezier() : GenericShape((ORDER_R + 1) * (ORDER_S + 1), 4) {}
61
62 IntegrationScheme create_ischeme() const {
63 IS_Q orig_ischeme;
64 orig_ischeme[0].set_order(2 * ORDER_R);
65 orig_ischeme[1].set_order(2 * ORDER_S);
66 // orig_ischeme[0].set_num(9);
67 // orig_ischeme[1].set_num(9);
68 IntegrationScheme ischeme;
69 for (int i = 0; i != orig_ischeme.get_num(); ++i) {
70 double xi[2];
71 double weight;
72 orig_ischeme.get_point(i, xi, weight);
73 xi[0] = .5 * (xi[0] + 1);
74 xi[1] = .5 * (xi[1] + 1);
75 weight *= .25;
76 ischeme.points.push_back(IntegrationScheme::Point(xi, weight));
77 }
78 return ischeme;
79 }
80
81 int get_edge_order(const int edge) const {
82 assert(edge >= 0 && edge < 4);
83 return (edge % 2 == 0 ? 2 * ORDER_R : 2 * ORDER_S); // XXX change
84 }
85
86 int get_edge_shape_order(const int edge) const {
87 assert(edge >= 0 && edge < 4);
88 return (edge % 2 == 0 ? ORDER_R : ORDER_S);
89 }
90
91 void get_xi_all_from_xi_edge(const int edge, const double xi, double xi_all[2]) const {
92 assert(edge >= 0 && edge < 4);
93 const double mapping[4][2] = {{xi, 0.}, {1., xi}, {1 - xi, 1.}, {0., 1 - xi}};
94 xi_all[0] = mapping[edge][0];
95 xi_all[1] = mapping[edge][1];
96 }
97
98 IntegrationScheme create_ischeme_edge(const int edge) const {
99 assert(edge >= 0 && edge < 4);
100 IS_L orig_ischeme;
101 orig_ischeme.set_order(get_edge_order(edge));
102
103 IntegrationScheme ischeme;
104 for (int i = 0; i != orig_ischeme.get_num(); ++i) {
105 double xi;
106 double weight;
107 orig_ischeme.get_point(i, xi, weight);
108 xi = .5 * (xi + 1);
109 weight *= 0.5;
110 const double mapping[4][2] = {{xi, 0.}, {1., xi}, {1 - xi, 1.}, {0., 1 - xi}};
111 ischeme.points.push_back(
112 IntegrationScheme::Point(mapping[edge][0], mapping[edge][1], weight));
113 }
114 return ischeme;
115 }
116
117 void get_edge_direction(const int edge, double edge_dir[2]) const {
118 assert(edge >= 0 && edge < 4);
119 switch (edge) {
120 case 0:
121 edge_dir[0] = 1;
122 edge_dir[1] = 0;
123 break;
124 case 1:
125 edge_dir[0] = 0;
126 edge_dir[1] = 1;
127 break;
128 case 2:
129 edge_dir[0] = -1;
130 edge_dir[1] = 0;
131 break;
132 case 3:
133 edge_dir[0] = 0;
134 edge_dir[1] = -1;
135 break;
136 }
137 }
138
139 void get_node_indices(const int sub_type, const int sub_index, b2linalg::Index& indices) const {
140 indices.clear();
141 if (sub_type == SUB_ALL) {
142 indices.reserve(num_nodes);
143 for (size_t i = 0; i != num_nodes; ++i) { indices.push_back(i); }
144 } else if (sub_type == SUB_EDGE_C0 || sub_type == SUB_EDGE_C1 || sub_type == SUB_EDGE_C2) {
145 const int edge = sub_index;
146 assert(edge >= 0 && edge < 4);
147 const int continuity =
148 (sub_type == SUB_EDGE_C0 ? 0 : (sub_type == SUB_EDGE_C1 ? 1 : 2));
149 if (edge % 2 == 0) {
150 const size_t n = num_nodes_edge_r * (1 + continuity);
151 indices.reserve(n);
152 for (size_t i = 0; i != n; ++i) {
153 const size_t ii = DATA::node_indices_edge_c2_r[edge / 2][i];
154 indices.push_back(ii);
155 }
156 } else {
157 const size_t n = num_nodes_edge_s * (1 + continuity);
158 indices.reserve(n);
159 for (size_t i = 0; i != n; ++i) {
160 const size_t ii = DATA::node_indices_edge_c2_s[edge / 2][i];
161 indices.push_back(ii);
162 }
163 }
164 } else if (
165 sub_type == SUB_VERTEX_C0 || sub_type == SUB_VERTEX_C1 || sub_type == SUB_VERTEX_C2) {
166 const int vertex = sub_index;
167 assert(vertex >= 0 && vertex < 4);
168 const int continuity =
169 (sub_type == SUB_VERTEX_C0 ? 0 : (sub_type == SUB_VERTEX_C1 ? 1 : 2));
170 const size_t n = (continuity == 0 ? 1 : (continuity == 1 ? 3 : 6));
171 indices.reserve(n);
172 for (size_t i = 0; i != n; ++i) {
173 const size_t ii = DATA::node_indices_vertex_c2[vertex][i];
174 indices.push_back(ii);
175 }
176 }
177 }
178
179 void compute_nodes_interpolation_straight_d0(
180 const double xi[2], b2linalg::Vector<double, b2linalg::Vdense>& N) const {
181 N.resize(4);
182 N[0] = (1 - xi[0]) * (1 - xi[1]);
183 N[1] = xi[0] * (1 - xi[1]);
184 N[2] = xi[0] * xi[1];
185 N[3] = (1 - xi[0]) * xi[1];
186 }
187
188 void compute_nodes_interpolation_straight_d1(
189 const double xi[2], b2linalg::Matrix<double, b2linalg::Mrectangle>& N) const {
190 N.resize(4, 3);
191
192 N[d00][0] = (1 - xi[0]) * (1 - xi[1]);
193 N[d00][1] = xi[0] * (1 - xi[1]);
194 N[d00][2] = xi[0] * xi[1];
195 N[d00][3] = (1 - xi[0]) * xi[1];
196
197 N[d10][0] = -(1 - xi[1]);
198 N[d10][1] = (1 - xi[1]);
199 N[d10][2] = xi[1];
200 N[d10][3] = -xi[1];
201
202 N[d01][0] = -(1 - xi[0]);
203 N[d01][1] = -xi[0];
204 N[d01][2] = xi[0];
205 N[d01][3] = (1 - xi[0]);
206 }
207
208 void compute_nodes_interpolation_straight_d2(
209 const double xi[2], b2linalg::Matrix<double, b2linalg::Mrectangle>& N) const {
210 N.resize(4, 6);
211
212 N[d00][0] = (1 - xi[0]) * (1 - xi[1]);
213 N[d00][1] = xi[0] * (1 - xi[1]);
214 N[d00][2] = xi[0] * xi[1];
215 N[d00][3] = (1 - xi[0]) * xi[1];
216
217 N[d10][0] = -(1 - xi[1]);
218 N[d10][1] = (1 - xi[1]);
219 N[d10][2] = xi[1];
220 N[d10][3] = -xi[1];
221
222 N[d01][0] = -(1 - xi[0]);
223 N[d01][1] = -xi[0];
224 N[d01][2] = xi[0];
225 N[d01][3] = (1 - xi[0]);
226
227 N[d20][0] = 0;
228 N[d20][1] = 0;
229 N[d20][2] = 0;
230 N[d20][3] = 0;
231
232 N[d11][0] = +1;
233 N[d11][1] = -1;
234 N[d11][2] = +1;
235 N[d11][3] = -1;
236
237 N[d02][0] = 0;
238 N[d02][1] = 0;
239 N[d02][2] = 0;
240 N[d02][3] = 0;
241 }
242};
243
244template <int ORDER_R, int ORDER_S, typename DATA>
245const size_t QuadBezier<ORDER_R, ORDER_S, DATA>::num_nodes_edge_c1_r;
246
247template <int ORDER_R, int ORDER_S, typename DATA>
248const size_t QuadBezier<ORDER_R, ORDER_S, DATA>::num_nodes_edge_c2_r;
249
250template <int ORDER_R, int ORDER_S, typename DATA>
251const size_t QuadBezier<ORDER_R, ORDER_S, DATA>::num_nodes_edge_c1_s;
252
253template <int ORDER_R, int ORDER_S, typename DATA>
254const size_t QuadBezier<ORDER_R, ORDER_S, DATA>::num_nodes_edge_c2_s;
255
256} // namespace b2000::klt
257
258#endif // B2ELEMENT_KLT_COMPUTE_SHAPE_Q_H_
b2tensor_calculus.H