b2api
B2000++ API Reference Manual, VERSION 4.6
 
Loading...
Searching...
No Matches
b2element_contact.H
1//------------------------------------------------------------------------
2// b2element_contact.H --
3//
4//
5// written by Mathias Doreille
6// Neda Ebrahimi Pour <neda.ebrahimipour@dlr.de>
7// Thomas Blome <thomas.blome@dlr.de>
8//
9// (c) 2009, 2016 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_CONTACT_H_
22#define B2ELEMENT_CONTACT_H_
23
24#include <limits>
25
26#include "b2element_continuum_integrate.H"
27#include "elements/properties/b2contact_mechanics.H"
28#include "model/b2domain.H"
29#include "model/b2element.H"
30#include "model/b2model.H"
31#include "model_imp/b2hnode.H"
32#include "utils/b2point_triangle_distance.H"
33#include "utils/b2tensor_calculus.H"
34
35namespace b2000 {
36
37typedef double darrar3[3];
38
39template <int N>
40struct Pow3 {
41 static const int value = Pow3<N - 1>::value * 3;
42};
43
44template <>
45struct Pow3<0> {
46 static const int value = 1;
47};
48
49template <int NB_LEVEL>
50class TriangleOnePointCompositIntegration {
51public:
52 static int nb_point[NB_LEVEL];
53 static darrar3 weight[NB_LEVEL][Pow3<NB_LEVEL>::value];
54
55 TriangleOnePointCompositIntegration() {
56 const darrar3 w1 = {1, 0, 0};
57 const darrar3 w2 = {0, 1, 0};
58 const darrar3 w3 = {0, 0, 1};
59 for (int i = 0; i < NB_LEVEL; ++i) {
60 const darrar3* ptr = init(i, w1, w2, w3, weight[i]);
61 nb_point[i] = ptr - weight[i];
62 }
63 }
64
65private:
66 static darrar3* init(
67 int level, const darrar3 w1, const darrar3 w2, const darrar3 w3, darrar3* ptr) {
68 if (level == 0) {
69 (*ptr)[0] = (w1[0] + w2[0] + w3[0]) / 3;
70 (*ptr)[1] = (w1[1] + w2[1] + w3[1]) / 3;
71 (*ptr)[2] = (w1[2] + w2[2] + w3[2]) / 3;
72 return ptr + 1;
73 } else {
74 const darrar3 w4 = {
75 (w1[0] + w2[0] + w3[0]) / 3, (w1[1] + w2[1] + w3[1]) / 3,
76 (w1[2] + w2[2] + w3[2]) / 3};
77 --level;
78 ptr = init(level, w1, w2, w4, ptr);
79 ptr = init(level, w2, w3, w4, ptr);
80 ptr = init(level, w3, w1, w4, ptr);
81 return ptr;
82 }
83 }
84 static TriangleOnePointCompositIntegration dumy;
85};
86
87template <int NB_LEVEL>
88TriangleOnePointCompositIntegration<NB_LEVEL> TriangleOnePointCompositIntegration<NB_LEVEL>::dumy;
89
90template <int NB_LEVEL>
91darrar3 TriangleOnePointCompositIntegration<NB_LEVEL>::weight[NB_LEVEL][Pow3<NB_LEVEL>::value];
92
93template <int NB_LEVEL>
94int TriangleOnePointCompositIntegration<NB_LEVEL>::nb_point[NB_LEVEL];
95
96// #define WITHOUT_CONSTRAINT 1
97
98// #ifndef WITHOUT_CONSTRAINT
99// #define USE_UNSYMMETRIC 1
100// #endif
101
102#define USE_UNSYMMETRIC 1
103// #define CHECK_DERIVATIVE 1
104
105class ElementContactT3 : public TypedElement<double> {
106public:
107 using node_type = node::HNode<
108 coordof::Coor<coordof::Trans<coordof::X, coordof::Y, coordof::Z>>,
109 coordof::Dof<coordof::DTrans<coordof::DX, coordof::DY, coordof::DZ>, coordof::Contact>>;
110
111 using type_t = ObjectTypeComplete<ElementContactT3, TypedElement<double>::type_t>;
112
113 ElementContactT3()
114#ifdef USE_UNSYMMETRIC
115 : TypedElement{false}
116#endif
117 {
118 }
119
120 void set_nodes(std::pair<int, Node* const*> nodes_) override {
121 if (nodes_.first != 3) { Exception() << "This element has 3 nodes." << THROW; }
122 for (int i = 0; i != 3; ++i) {
123 nodes[i] = nodes_.second[i];
124 if (!node_type::is_dof_subset_s(nodes[i])) {
125 Exception() << "This element cannot accept node of type " << typeid(*nodes[i])
126 << ". The only accepted node types are " << typeid(node_type) << "."
127 << THROW;
128 }
129 }
130 }
131
132 std::pair<int, Node* const*> get_nodes() const override {
133 return std::pair<int, Node* const*>(3, nodes);
134 }
135
136 void set_property(ElementProperty* property_) override {
137 /*property = dynamic_cast<ContactMechanics*>(property_);
138 if (property == 0)
139 TypeError() << THROW;
140 */
141 }
142
143 const ElementProperty* get_property() const override { return property; }
144
145 void init(Model& model) override {}
146
147 void get_dof_numbering(b2linalg::Index& dof_numbering) override {
148 dof_numbering.resize(12);
149 size_t* ptr = &dof_numbering[0];
150 for (int k = 0; k != 3; ++k) { ptr = node_type::get_global_dof_numbering_s(ptr, nodes[k]); }
151 }
152
153 const std::vector<VariableInfo> get_value_info() const override {
154 return std::vector<VariableInfo>(2, Element::nonlinear);
155 }
156
157#ifdef OLD_IMPLE
158 void get_value(
159 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
160 const double time, const double delta_time,
161 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof,
162 GradientContainer* gradient_container, SolverHints* solver_hints,
163 b2linalg::Vector<double, b2linalg::Vdense>& value,
164 std::vector<b2linalg::Matrix<
165 double,
166#ifdef USE_UNSYMMETRIC
167 b2linalg::Mrectangle
168#else
169 b2linalg::Mpacked
170#endif
171 >>& d_value_d_dof,
172 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_time) {
173
174 typedef b2linalg::Matrix<
175 double,
176#ifdef USE_UNSYMMETRIC
177 b2linalg::Mrectangle
178#else
179 b2linalg::Mpacked
180#endif
181 >
182 matrix_t;
183
184 const double epsilon1 = 1e1;
185 const double epsilon2 = 1e2;
186
187 // The coordinate at each nodes
188 double R[3][3];
189 // The displacement at each nodes
190 double u[3][3];
191 // The pressure at each nodes
192 double lambda[3];
193
194 if (dof.is_null()) {
195 for (int k = 0; k != 3; ++k) { node_type::get_coor_s(R[k], nodes[k]); }
196 std::fill_n(u[0], 9, 0);
197 std::fill_n(lambda, 3, 0);
198 } else {
199 const double* ptr = &dof(0, 0);
200 for (int k = 0; k != 3; ++k) {
201 node_type::get_coor_s(R[k], nodes[k]);
202 u[k][0] = *ptr++;
203 u[k][1] = *ptr++;
204 u[k][2] = *ptr++;
205 lambda[k] = *ptr++;
206 }
207 }
208
209 double surface_initial;
210 // double normal_scale;
211 {
212 double J[2][3];
213 for (int i = 0; i != 3; ++i) {
214 J[0][i] = -R[0][i] + R[1][i];
215 J[1][i] = -R[0][i] + R[2][i];
216 }
217 double normal[3];
218 outer_product_3(J[0], J[1], normal);
219 surface_initial = 0.5 * std::sqrt(dot_3(normal, normal));
220 /*
221 for (int i = 0; i != 3; ++i) {
222 J[0][i] += -u[0][i] + u[1][i];
223 J[1][i] += -u[0][i] + u[2][i];
224 }
225 outer_product_3(J[0], J[1], normal);
226 normal_scale = 1 / normalise_3(normal);
227 */
228 }
229
230 if (!value.is_null()) {
231 value.resize(12);
232 value.set_zero();
233 }
234
235 bool comput_d_value_d_dof = !d_value_d_dof.empty() && !d_value_d_dof[0].is_null();
236
237 matrix_t matrix_tmp;
238 matrix_t& matrix =
239 0 /*(!comput_d_value_d_dof && !value.is_null() && dof.size2() > 1)*/ ? matrix_tmp
240 : d_value_d_dof
241 [0];
242 if (&matrix == &matrix_tmp) { comput_d_value_d_dof = true; }
243
244 if (comput_d_value_d_dof) {
245#ifdef USE_UNSYMMETRIC
246 matrix.resize(12, 12);
247#else
248 matrix.resize(12, 12, true);
249#endif
250 matrix.set_zero();
251 /*
252 // d_normal_d_u[k][i][j] = d Ni / d u[k][j]
253 double d_normal_d_u[3][3][3];
254 d_normal_d_u[0][0][0] = 0;
255 d_normal_d_u[0][0][1] = normal_scale * (u[1][2] - u[2][2]);
256 d_normal_d_u[0][0][2] = normal_scale * (u[2][1] - u[1][1]);
257 d_normal_d_u[1][0][0] = 0;
258 d_normal_d_u[1][0][1] = normal_scale * (u[2][2] - u[0][2]);
259 d_normal_d_u[1][0][2] = normal_scale * (u[0][1] - u[2][1]);
260 d_normal_d_u[2][0][0] = 0;
261 d_normal_d_u[2][0][1] = normal_scale * (u[0][2] - u[1][2]);
262 d_normal_d_u[2][0][2] = normal_scale * (u[1][1] - u[0][1]);
263
264 d_normal_d_u[0][1][0] = normal_scale * (u[2][2] - u[1][2]);
265 d_normal_d_u[0][1][1] = 0;
266 d_normal_d_u[0][1][2] = normal_scale * (u[1][0] - u[2][0]);
267 d_normal_d_u[1][1][0] = normal_scale * (u[0][2] - u[2][2]);
268 d_normal_d_u[1][1][1] = 0;
269 d_normal_d_u[1][1][2] = normal_scale * (u[2][0] - u[0][0]);
270 d_normal_d_u[2][1][0] = normal_scale * (u[1][2] - u[0][2]);
271 d_normal_d_u[2][1][1] = 0;
272 d_normal_d_u[2][1][2] = normal_scale * (u[0][0] - u[1][0]);
273
274 d_normal_d_u[0][2][0] = normal_scale * (u[1][1] - u[2][1]);
275 d_normal_d_u[0][2][1] = normal_scale * (u[2][0] - u[1][0]);
276 d_normal_d_u[0][2][2] = 0;
277 d_normal_d_u[1][2][0] = normal_scale * (u[2][1] - u[0][1]);
278 d_normal_d_u[1][2][1] = normal_scale * (u[0][0] - u[2][0]);
279 d_normal_d_u[1][2][2] = 0;
280 d_normal_d_u[2][2][0] = normal_scale * (u[0][1] - u[1][1]);
281 d_normal_d_u[2][2][1] = normal_scale * (u[1][0] - u[0][0]);
282 d_normal_d_u[2][2][2] = 0;
283 */
284 }
285
286 IS_T integration;
287 integration.set_num(3);
288 for (int p = 0; p != integration.get_num(); ++p) {
289 IP_ID ip_id;
290 double weight;
291 integration.get_point(p, ip_id, weight);
292 weight *= surface_initial;
293 const double shape[3] = {1 - ip_id[0] - ip_id[1], ip_id[0], ip_id[1]};
294 double x[3] = {0, 0, 0};
295 double l = 0;
296 for (int k = 0; k != 3; ++k) {
297 l += shape[k] * lambda[k];
298 for (int i = 0; i != 3; ++i) { x[i] += shape[k] * (u[k][i] + R[k][i]); }
299 }
300 const double g = x[2]; // the contact surface is z > 0
301 const double gnormal[3] = {0, 0, -1}; // The normal of the contact
302 double d_gnormal_d_u[3][3][3];
303 std::fill_n(d_gnormal_d_u[0][0], 27, 0);
304 const double d_g_d_u[3][3] = {
305 {
306 0,
307 0,
308 0,
309 },
310 {0, 0, 0},
311 { shape[0],
312 shape[1],
313 shape[2] }};
314 const double w = epsilon1 * ((g + l) / 2 - std::sqrt((g - l) * (g - l) / 4 + epsilon2));
315 // std::cout << "p=" << p << ", x=" << x[0] << " " << x[1] << " " << x[2] << ", g=" << g
316 // << ", l=" << l << ", w=" << w << std::endl;
317#ifdef WITHOUT_CONSTRAINT
318 if (g < 1e-4 && l > 1e-4) {
319 if (!value.is_null()) {
320 double* ptr = &value[0];
321 for (int k = 0; k != 3; ++k, ++ptr) {
322 for (int i = 0; i != 3; ++i, ++ptr) {
323 for (int j = 0; j != 3; ++j) {
324 *ptr -= weight * shape[j] * d_g_d_u[i][k] * lambda[j];
325 }
326 }
327 *ptr -= weight * shape[k] * g;
328 }
329 }
330 if (comput_d_value_d_dof) {
331 double* ptr = &d_value_d_dof[0](0, 0);
332 for (int ki = 0; ki != 3; ++ki) {
333 for (int i = 0; i != 3; ++i) {
334 for (int kj = 0; kj <= ki; ++kj) {
335 if (ki == kj) {
336 for (int j = 0; j <= i; ++j) { ptr++; }
337 } else {
338 for (int j = 0; j != 3; ++j) { ptr++; }
339 *ptr++ -= weight * shape[kj] * d_g_d_u[i][ki];
340 }
341 }
342 }
343 for (int kj = 0; kj <= ki; ++kj) {
344 for (int j = 0; j != 3; ++j) {
345 *ptr++ -= weight * shape[ki] * d_g_d_u[j][kj];
346 }
347 ptr++;
348 }
349 }
350 }
351 } else {
352 if (!value.is_null()) {
353 for (int ki = 0; ki != 3; ++ki) {
354 value[ki * 4 + 3] += 1e-8 * weight * shape[ki] * l;
355 }
356 }
357 if (comput_d_value_d_dof) {
358 for (int ki = 0; ki != 3; ++ki) {
359 d_value_d_dof[0](ki * 4 + 3, ki * 4 + 3) +=
360 1e-8 * weight * shape[ki] * shape[ki];
361 }
362 }
363 }
364#else
365 if (!value.is_null()) {
366 double* ptr = &value[0];
367 for (int k = 0; k != 3; ++k, ++ptr) {
368 for (int i = 0; i != 3; ++i, ++ptr) {
369 *ptr += weight * gnormal[i] * shape[k] * l;
370 }
371 *ptr += weight * shape[k] * w;
372 }
373 }
374 if (comput_d_value_d_dof) {
375 const double d_w_d_g =
376 epsilon1 * 0.5 * (1 - (g - l) / std::sqrt((g - l) * (g - l) + 4 * epsilon2));
377 const double d_w_d_l =
378 epsilon1 * 0.5 * (1 + (g - l) / std::sqrt((g - l) * (g - l) + 4 * epsilon2));
379 double* ptr = &matrix(0, 0);
380#ifdef USE_UNSYMMETRIC
381 for (int kj = 0; kj != 3; ++kj) {
382 for (int j = 0; j != 3; ++j) {
383 for (int ki = 0; ki != 3; ++ki) {
384 for (int i = 0; i != 3; ++i) {
385 *ptr++ += weight * l * d_gnormal_d_u[kj][i][j] * shape[ki];
386 }
387 *ptr++ += weight * shape[ki] * d_w_d_g * d_g_d_u[j][kj];
388 }
389 }
390 for (int ki = 0; ki != 3; ++ki) {
391 for (int i = 0; i != 3; ++i) {
392 *ptr++ += weight * shape[kj] * gnormal[i] * shape[ki];
393 }
394 *ptr++ += weight * shape[ki] * d_w_d_l * shape[kj];
395 }
396 }
397 /*
398 for (int i = 0; i != 12; ++i)
399 for (int j = i + 1; j < 12; ++j)
400 if (d_value_d_dof[0](i, j) != d_value_d_dof[0](j, i))
401 std::cout << get_object_name() << " " << i << " " << j << " " <<
402 d_value_d_dof[0](i, j) << " " << d_value_d_dof[0](j, i) << std::endl;
403 */
404#else
405 for (int ki = 0; ki != 3; ++ki) {
406 for (int i = 0; i != 3; ++i) {
407 for (int kj = 0; kj <= ki; ++kj) {
408 const int j_max = kj == ki ? i : 2;
409 for (int j = 0; j <= j_max; ++j) {
410 *ptr++ += weight * l * 0.5
411 * (d_gnormal_d_u[kj][i][j] * shape[ki]
412 + d_gnormal_d_u[ki][j][i] * shape[kj]);
413 }
414 if (kj != ki) {
415 *ptr++ += weight * shape[kj] * 0.5
416 * (gnormal[i] * shape[ki] + d_w_d_g * d_g_d_u[i][ki]);
417 }
418 }
419 }
420 for (int kj = 0; kj <= ki; ++kj) {
421 for (int j = 0; j != 3; ++j) {
422 *ptr++ += weight * shape[ki] * 0.5
423 * (d_w_d_g * d_g_d_u[j][kj] + gnormal[j] * shape[kj]);
424 }
425 *ptr++ += weight * shape[ki] * d_w_d_l * shape[kj];
426 }
427 }
428#endif
429 }
430#endif
431 }
432
433 // if (!value.is_null())
434 // for (int i = 1; i < dof.size2(); ++i)
435 // value += matrix * dof[i];
436 for (int i = 1; i < d_value_d_dof.size(); ++i) {
437 if (!d_value_d_dof[i].is_null()) {
438#ifdef USE_UNSYMMETRIC
439 d_value_d_dof[i].resize(12, 12);
440#else
441 d_value_d_dof[i].resize(12, 12, true);
442#endif
443 d_value_d_dof[i].set_zero(); // = matrix;
444 }
445 }
446 }
447#endif
448
449 struct ElemCoor {
450 double coor[3][3];
451 };
452
453 void get_axe_aligned_bounding_box(
454 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof, int dim, double* min,
455 double* max) override {
456 for (size_t i = 0; i != 3; ++i) {
457 min[i] = std::numeric_limits<double>::max();
458 max[i] = std::numeric_limits<double>::lowest();
459 }
460 for (size_t p = 0; p != 3; ++p) {
461 double coor[3];
462 node_type::get_coor_s(coor, nodes[p]);
463 size_t dn[4];
464 node_type::get_global_dof_numbering_s(dn, nodes[p]);
465 assert(dim <= 3);
466 for (int i = 0; i != dim; ++i) {
467 coor[i] += dof(dn[i], 0);
468 min[i] = std::min(min[i], coor[i]);
469 max[i] = std::max(max[i], coor[i]);
470 }
471 }
472 }
473
474 void get_value(
475 Model& model, const bool linear, const EquilibriumSolution equilibrium_solution,
476 const double time, const double delta_time,
477 const b2linalg::Matrix<double, b2linalg::Mrectangle_constref>& dof,
478 GradientContainer* gradient_container, SolverHints* solver_hints,
479 b2linalg::Index& dof_numbering, b2linalg::Vector<double, b2linalg::Vdense>& value,
480 const std::vector<bool>& d_value_d_dof_flags,
481 std::vector<b2linalg::Matrix<
482 double,
483#ifdef USE_UNSYMMETRIC
484 b2linalg::Mrectangle
485#else
486 b2linalg::Mpacked
487#endif
488 >>& d_value_d_dof,
489 b2linalg::Vector<double, b2linalg::Vdense>& d_value_d_time) override {
490 const double epsilon = 1e-0;
491 const double dist_box = 1;
492 const int integration_level = 1;
493 std::vector<Element*> element_list;
494 model.get_domain().get_elements_of_same_type_near(this, dist_box, element_list);
495
496 double coor[3][3];
497 double normal[3];
498
499 if (element_list.empty()) {
500 get_dof_numbering(dof_numbering);
501 for (int p = 0; p != 3; ++p) { node_type::get_coor_s(coor[p], nodes[p]); }
502 double normal[3];
503 {
504 double a[3] = {
505 coor[1][0] - coor[0][0], coor[1][1] - coor[0][1], coor[1][2] - coor[0][2]};
506 double b[3] = {
507 coor[2][0] - coor[0][0], coor[2][1] - coor[0][1], coor[2][2] - coor[0][2]};
508 outer_product_3(a, b, normal);
509 }
510 const double area =
511 normalise_3(normal) / (6 * composit_integration.nb_point[integration_level]);
512
513 for (int i = 0; i != 3; ++i) { dof_numbering[i] = dof_numbering[i * 4 + 3]; }
514 dof_numbering.resize(3);
515 if (!value.is_null()) {
516 value.resize(3);
517 for (int i = 0; i != 3; ++i) { value[i] = area * dof(dof_numbering[i], 0); }
518 }
519 if (!d_value_d_dof_flags.empty() && d_value_d_dof_flags[0]) {
520 d_value_d_dof[0].resize(
521 3, 3
522#ifndef USE_UNSYMMETRIC
523 ,
524 true
525#endif
526 );
527 d_value_d_dof[0].set_zero();
528 for (int i = 0; i != 3; ++i) { d_value_d_dof[0](i, i) = area; }
529 }
530 for (size_t i = 1; i < d_value_d_dof_flags.size(); ++i) {
531 if (d_value_d_dof_flags[i]) {
532 d_value_d_dof[i].resize(
533 dof_numbering.size(), dof_numbering.size()
534#ifndef USE_UNSYMMETRIC
535 ,
536 true
537#endif
538 );
539 d_value_d_dof[i].set_zero();
540 }
541 }
542 if (!d_value_d_time.is_null()) {
543 d_value_d_time.resize(dof_numbering.size());
544 d_value_d_time.set_zero();
545 }
546 return;
547 }
548
549 std::vector<ElemCoor> element_list_coor(element_list.size());
550 for (size_t e = 0; e != element_list_coor.size(); ++e) {
551 element_list[e]->get_dof_numbering(dof_numbering);
552 Node* const* n = element_list[e]->get_nodes().second;
553 for (int p = 0; p != 3; ++p) {
554 node_type::get_coor_s(element_list_coor[e].coor[p], n[p]);
555 for (int i = 0; i != 3; ++i) {
556 element_list_coor[e].coor[p][i] += dof(dof_numbering[p * 4 + i], 0);
557 }
558 }
559 }
560
561 get_dof_numbering(dof_numbering);
562
563 for (int p = 0; p != 3; ++p) { node_type::get_coor_s(coor[p], nodes[p]); }
564 {
565 double a[3] = {
566 coor[1][0] - coor[0][0], coor[1][1] - coor[0][1], coor[1][2] - coor[0][2]};
567 double b[3] = {
568 coor[2][0] - coor[0][0], coor[2][1] - coor[0][1], coor[2][2] - coor[0][2]};
569 outer_product_3(a, b, normal);
570 }
571 const double area =
572 normalise_3(normal) / (2 * composit_integration.nb_point[integration_level]);
573
574 for (int p = 0; p != 3; ++p) {
575 for (int i = 0; i != 3; ++i) { coor[p][i] += dof(dof_numbering[p * 4 + i], 0); }
576 }
577 {
578 double a[3] = {
579 coor[1][0] - coor[0][0], coor[1][1] - coor[0][1], coor[1][2] - coor[0][2]};
580 double b[3] = {
581 coor[2][0] - coor[0][0], coor[2][1] - coor[0][1], coor[2][2] - coor[0][2]};
582 outer_product_3(a, b, normal);
583 }
584 const double inv_norme_normal = 1 / normalise_3(normal);
585
586 // initialise the output vector and matrix
587 if (!value.is_null()) {
588 value.resize(12);
589 value.set_zero();
590 }
591
592 // d_normal_d_u[k][i][j] = d Nk / d u[i][j]
593 double d_normal_d_u[3][3][3];
594 if (!d_value_d_dof_flags.empty() && d_value_d_dof_flags[0]) {
595 d_value_d_dof[0].resize(
596 12, 12
597#ifndef USE_UNSYMMETRIC
598 ,
599 true
600#endif
601 );
602 d_value_d_dof[0].set_zero();
603
604 double d_n_d_u[3][3][3];
605 d_n_d_u[0][0][0] = 0;
606 d_n_d_u[0][0][1] = coor[1][2] - coor[2][2];
607 d_n_d_u[0][0][2] = coor[2][1] - coor[1][1];
608 d_n_d_u[0][1][0] = 0;
609 d_n_d_u[0][1][1] = coor[2][2] - coor[0][2];
610 d_n_d_u[0][1][2] = coor[0][1] - coor[2][1];
611 d_n_d_u[0][2][0] = 0;
612 d_n_d_u[0][2][1] = coor[0][2] - coor[1][2];
613 d_n_d_u[0][2][2] = coor[1][1] - coor[0][1];
614
615 d_n_d_u[1][0][0] = coor[2][2] - coor[1][2];
616 d_n_d_u[1][0][1] = 0;
617 d_n_d_u[1][0][2] = coor[1][0] - coor[2][0];
618 d_n_d_u[1][1][0] = coor[0][2] - coor[2][2];
619 d_n_d_u[1][1][1] = 0;
620 d_n_d_u[1][1][2] = coor[2][0] - coor[0][0];
621 d_n_d_u[1][2][0] = coor[1][2] - coor[0][2];
622 d_n_d_u[1][2][1] = 0;
623 d_n_d_u[1][2][2] = coor[0][0] - coor[1][0];
624
625 d_n_d_u[2][0][0] = coor[1][1] - coor[2][1];
626 d_n_d_u[2][0][1] = coor[2][0] - coor[1][0];
627 d_n_d_u[2][0][2] = 0;
628 d_n_d_u[2][1][0] = coor[2][1] - coor[0][1];
629 d_n_d_u[2][1][1] = coor[0][0] - coor[2][0];
630 d_n_d_u[2][1][2] = 0;
631 d_n_d_u[2][2][0] = coor[0][1] - coor[1][1];
632 d_n_d_u[2][2][1] = coor[1][0] - coor[0][0];
633 d_n_d_u[2][2][2] = 0;
634
635 for (int i = 0; i != 3; ++i) {
636 for (int j = 0; j != 3; ++j) {
637 double tmp = 0;
638 for (int l = 0; l != 3; ++l) { tmp += normal[l] * d_n_d_u[l][j][i]; }
639 for (int k = 0; k != 3; ++k) {
640 d_normal_d_u[k][j][i] =
641 inv_norme_normal * (d_n_d_u[k][j][i] - normal[k] * tmp);
642 }
643 }
644 }
645
646#ifdef CHECK_DERIVATIVE
647 {
648 const long double h = 0.0000001;
649 for (int j = 0; j != 3; ++j) {
650 for (int i = 0; i != 3; ++i) {
651 long double coor_tmp[3][3];
652 std::copy(coor[0], coor[3], coor_tmp[0]);
653 coor_tmp[j][i] += h;
654 long double a_tmp[3] = {
655 coor_tmp[1][0] - coor_tmp[0][0], coor_tmp[1][1] - coor_tmp[0][1],
656 coor_tmp[1][2] - coor_tmp[0][2]};
657 long double b_tmp[3] = {
658 coor_tmp[2][0] - coor_tmp[0][0], coor_tmp[2][1] - coor_tmp[0][1],
659 coor_tmp[2][2] - coor_tmp[0][2]};
660 long double normal_tmp[3];
661 outer_product_3(a_tmp, b_tmp, normal_tmp);
662 normalise_3(normal_tmp);
663 for (int k = 0; k != 3; ++k) {
664 long double v = (normal_tmp[k] - normal[k]) / h;
665 if (std::abs(d_normal_d_u[k][j][i] - v) > 1e-5) {
666 std::cout << "d_normal_d_u[" << k << "][" << j << "][" << i
667 << "]=" << d_normal_d_u[k][j][i] << ", n=" << v
668 << std::endl;
669 }
670 }
671 }
672 }
673 }
674#endif
675 }
676
677 std::map<size_t, size_t> list_contact_element;
678 for (int k = 0; k != composit_integration.nb_point[integration_level]; ++k) {
679 const double* weight = composit_integration.weight[integration_level][k];
680 double point[3] = {};
681 for (int p = 0; p != 3; ++p) {
682 for (int i = 0; i != 3; ++i) { point[i] += weight[p] * coor[p][i]; }
683 }
684 double dist = 1e100;
685 double r, s;
686 size_t elem = 0;
687 for (size_t i = 0; i != element_list.size(); ++i) {
688 double r_tmp, s_tmp;
689 const double dist_tmp = point_triangle_distance3d(
690 point, element_list_coor[i].coor[0], element_list_coor[i].coor[1],
691 element_list_coor[i].coor[2], r_tmp, s_tmp);
692 if (dist_tmp < dist) {
693 dist = dist_tmp;
694 r = r_tmp;
695 s = s_tmp;
696 elem = i;
697 }
698 }
699 if (dist > dist_box * dist_box) { continue; }
700
701 size_t dof_numbering_idx = list_contact_element[elem];
702 if (dof_numbering_idx == 0) {
703 list_contact_element[elem] = dof_numbering_idx = dof_numbering.size();
704 b2linalg::Index dof_numbering_tmp;
705 element_list[elem]->get_dof_numbering(dof_numbering_tmp);
706 for (int i = 0; i != 3; ++i) {
707 dof_numbering.insert(
708 dof_numbering.end(), dof_numbering_tmp.begin() + i * 4,
709 dof_numbering_tmp.begin() + i * 4 + 3);
710 }
711 if (!value.is_null()) { value.resize(dof_numbering.size()); }
712 if (!d_value_d_dof_flags.empty() && d_value_d_dof_flags[0]) {
713 d_value_d_dof[0].resize(
714 dof_numbering.size(), dof_numbering.size()
715#ifndef USE_UNSYMMETRIC
716 ,
717 true
718#endif
719 );
720 }
721 }
722
723 double lambda = 0;
724 double g = 0;
725 const double c_weight[3] = {1 - r - s, r, s};
726 for (int p = 0; p != 3; ++p) {
727 lambda += weight[p] * dof(dof_numbering[p * 4 + 3], 0);
728 g += (element_list_coor[elem].coor[0][p] * c_weight[0]
729 + element_list_coor[elem].coor[1][p] * r
730 + element_list_coor[elem].coor[2][p] * s - point[p])
731 * normal[p];
732 }
733 const double w =
734 area * ((g + lambda) / 2 - std::sqrt((g - lambda) * (g - lambda) / 4 + epsilon));
735 const double lambda_area = area * lambda;
736
737 /* {
738 std::cout << "point=" << point[0] << ", " << point[1] << ", " << point[2] <<
739 std::endl; std::cout << "normal=" << normal[0] << ", " << normal[1] << ", " << normal[2]
740 << std::endl; std::cout << "g=" << g << ", lambda=" << lambda << ", w=" << w / area <<
741 std::endl;
742 } */
743
744 if (!value.is_null()) {
745 for (int p = 0; p != 3; ++p) {
746 for (int i = 0; i != 3; ++i) {
747 const double v = lambda_area * normal[i];
748 value[p * 4 + i] += v * weight[p];
749 value[dof_numbering_idx + p * 3 + i] -= v * c_weight[p];
750 }
751 value[p * 4 + 3] -= weight[p] * w;
752 }
753 }
754
755 if (!d_value_d_dof_flags.empty() && d_value_d_dof_flags[0]) {
756 double d_r[4][3];
757 double d_s[4][3];
758 point_triangle_distance3d(
759 point, element_list_coor[elem].coor[0], element_list_coor[elem].coor[1],
760 element_list_coor[elem].coor[2], r, s, d_r, d_s);
761
762#ifdef CHECK_DERIVATIVE
763 {
764 const double h = 0.000001;
765 for (int i = 0; i != 3; ++i) {
766 double point_tmp[3] = {point[0], point[1], point[2]};
767 point_tmp[i] += h;
768 double r1, s1;
769 point_triangle_distance3d(
770 point_tmp, element_list_coor[elem].coor[0],
771 element_list_coor[elem].coor[1], element_list_coor[elem].coor[2], r1,
772 s1);
773 double rtmp = (r1 - r) / h;
774 double stmp = (s1 - s) / h;
775 if (std::abs(d_r[0][i] - rtmp) > 1e-5) {
776 std::cout << "d_r[0][" << i << "]=" << d_r[0][i] << ", n=" << rtmp
777 << std::endl;
778 }
779 if (std::abs(d_s[0][i] - stmp) > 1e-5) {
780 std::cout << "d_s[0][" << i << "]=" << d_s[0][i] << ", n=" << stmp
781 << std::endl;
782 }
783 }
784 for (int i = 0; i != 3; ++i) {
785 for (int j = 0; j != 3; ++j) {
786 double coor_tmp[3][3];
787 for (int k = 0; k != 3; ++k) {
788 for (int l = 0; l != 3; ++l) {
789 coor_tmp[k][l] = element_list_coor[elem].coor[k][l];
790 }
791 }
792 coor_tmp[i][j] += h;
793 double r1, s1;
794 point_triangle_distance3d(
795 point, coor_tmp[0], coor_tmp[1], coor_tmp[2], r1, s1);
796 double rtmp = (r1 - r) / h;
797 double stmp = (s1 - s) / h;
798 if (std::abs(rtmp - d_r[i + 1][j]) > 1e-5) {
799 std::cout << "d_r[" << i + 1 << "][" << j << "]=" << d_r[i + 1][j]
800 << ", n=" << rtmp << std::endl;
801 }
802 if (std::abs(stmp - d_s[i + 1][j]) > 1e-5) {
803 std::cout << "d_s[" << i + 1 << "][" << j << "]=" << d_s[i + 1][j]
804 << ", n=" << stmp << std::endl;
805 }
806 }
807 }
808 }
809#endif
810
811 double d_c_weight[3][6][3];
812 for (int p = 0; p != 3; ++p) {
813 for (int i = 0; i != 3; ++i) {
814 d_c_weight[0][p][i] = -weight[p] * (d_r[0][i] + d_s[0][i]);
815 d_c_weight[1][p][i] = weight[p] * d_r[0][i];
816 d_c_weight[2][p][i] = weight[p] * d_s[0][i];
817 d_c_weight[0][p + 3][i] = -d_r[p + 1][i] - d_s[p + 1][i];
818 d_c_weight[1][p + 3][i] = d_r[p + 1][i];
819 d_c_weight[2][p + 3][i] = d_s[p + 1][i];
820 }
821 }
822
823 const double tmp1 = std::sqrt((lambda - g) * (lambda - g) + 4 * epsilon);
824 const double d_w_d_g = area * (lambda - g + tmp1) / (2 * tmp1);
825 const double d_w_d_lambda = area * (g - lambda + tmp1) / (2 * tmp1);
826 // std::cout << "d_w_d_g=" << d_w_d_g / area << ", d_w_d_lambda=" << d_w_d_lambda /
827 // area << std::endl;
828#ifdef CHECK_DERIVATIVE
829 {
830 const double h = 0.000001;
831 double g_tmp = g + h;
832 double w_tmp =
833 area
834 * ((g_tmp + lambda) / 2
835 - std::sqrt((g_tmp - lambda) * (g_tmp - lambda) / 4 + epsilon));
836 double d_w_d_g_tmp = (w_tmp - w) / h;
837 if (std::abs(d_w_d_g_tmp - d_w_d_g) > 1e-5) {
838 std::cout << "d_w_d_g=" << d_w_d_g << ", n=" << d_w_d_g_tmp << std::endl;
839 }
840 double lambda_tmp = lambda + h;
841 w_tmp = area
842 * ((g + lambda_tmp) / 2
843 - std::sqrt((g - lambda_tmp) * (g - lambda_tmp) / 4 + epsilon));
844 double d_w_d_lambda_tmp = (w_tmp - w) / h;
845 if (std::abs(d_w_d_lambda_tmp - d_w_d_lambda) > 1e-5) {
846 std::cout << "d_w_d_lambda=" << d_w_d_lambda << ", n=" << d_w_d_lambda_tmp
847 << std::endl;
848 }
849 }
850#endif
851 double d_g_d_u[6][3] = {};
852 for (int i = 0; i != 3; ++i) {
853 const double tmp = element_list_coor[elem].coor[0][i] * c_weight[0]
854 + element_list_coor[elem].coor[1][i] * r
855 + element_list_coor[elem].coor[2][i] * s - point[i];
856 for (int j = 0; j != 3; ++j) {
857 for (int k = 0; k != 3; ++k) {
858 d_g_d_u[j][k] += tmp * d_normal_d_u[i][j][k];
859 }
860 d_g_d_u[j][i] -= normal[i] * weight[j];
861 d_g_d_u[j + 3][i] += normal[i] * c_weight[j];
862 }
863 const double tmp1 = element_list_coor[elem].coor[i][0] * normal[0]
864 + element_list_coor[elem].coor[i][1] * normal[1]
865 + element_list_coor[elem].coor[i][2] * normal[2];
866 for (int j = 0; j != 6; ++j) {
867 for (int k = 0; k != 3; ++k) {
868 d_g_d_u[j][k] += tmp1 * d_c_weight[i][j][k];
869 }
870 }
871 }
872
873#ifdef CHECK_DERIVATIVE
874 {
875 const double h = 0.0000001;
876 for (int j = 0; j != 3; ++j) {
877 for (int i = 0; i != 3; ++i) {
878 double coor_tmp[3][3];
879 for (int k = 0; k != 3; ++k) {
880 for (int l = 0; l != 3; ++l) { coor_tmp[k][l] = coor[k][l]; }
881 }
882 coor_tmp[j][i] += h;
883 double point_tmp[3] = {};
884 for (int k = 0; k != 3; ++k) {
885 for (int l = 0; l != 3; ++l) {
886 point_tmp[l] += weight[k] * coor_tmp[k][l];
887 }
888 }
889 double normal_tmp[3];
890 {
891 double a[3] = {
892 coor_tmp[1][0] - coor_tmp[0][0],
893 coor_tmp[1][1] - coor_tmp[0][1],
894 coor_tmp[1][2] - coor_tmp[0][2]};
895 double b[3] = {
896 coor_tmp[2][0] - coor_tmp[0][0],
897 coor_tmp[2][1] - coor_tmp[0][1],
898 coor_tmp[2][2] - coor_tmp[0][2]};
899 outer_product_3(a, b, normal_tmp);
900 }
901 normalise_3(normal_tmp);
902
903 double rtmp, stmp;
904 point_triangle_distance3d(
905 point_tmp, element_list_coor[elem].coor[0],
906 element_list_coor[elem].coor[1], element_list_coor[elem].coor[2],
907 rtmp, stmp);
908 double gtmp = 0;
909 const double c_weight_tmp[3] = {1 - rtmp - stmp, rtmp, stmp};
910 for (int p = 0; p != 3; ++p) {
911 gtmp += (element_list_coor[elem].coor[0][p] * c_weight_tmp[0]
912 + element_list_coor[elem].coor[1][p] * rtmp
913 + element_list_coor[elem].coor[2][p] * stmp - point_tmp[p])
914 * normal_tmp[p];
915 }
916 double g1 = (gtmp - g) / h;
917 if (std::abs(g1 - d_g_d_u[j][i]) > 1e-5) {
918 std::cout << "d_g_d_u[" << j << "][" << i << "]=" << d_g_d_u[j][i]
919 << ", n=" << g1 << std::endl;
920 }
921 }
922 }
923 for (int j = 0; j != 3; ++j) {
924 for (int i = 0; i != 3; ++i) {
925 double coor_tmp[3][3];
926 for (int k = 0; k != 3; ++k) {
927 for (int l = 0; l != 3; ++l) {
928 coor_tmp[k][l] = element_list_coor[elem].coor[k][l];
929 }
930 }
931 coor_tmp[j][i] += h;
932 double rtmp, stmp;
933 point_triangle_distance3d(
934 point, coor_tmp[0], coor_tmp[1], coor_tmp[2], rtmp, stmp);
935 double gtmp = 0;
936 const double c_weight[3] = {1 - rtmp - stmp, rtmp, stmp};
937 for (int p = 0; p != 3; ++p) {
938 gtmp += (coor_tmp[0][p] * c_weight[0] + coor_tmp[1][p] * rtmp
939 + coor_tmp[2][p] * stmp - point[p])
940 * normal[p];
941 }
942 double g1 = (gtmp - g) / h;
943 if (std::abs(g1 - d_g_d_u[j + 3][i]) > 1e-5) {
944 std::cout << "d_g_d_u[" << j + 3 << "][" << i
945 << "]=" << d_g_d_u[j + 3][i] << ", n=" << g1 << std::endl;
946 }
947 }
948 }
949 }
950#endif
951
952 for (int pi = 0; pi != 3; ++pi) {
953 for (int i = 0; i != 3; ++i) {
954 const int idx = pi * 4 + i;
955 for (int pj = 0; pj != 3; ++pj) {
956 for (int j = 0; j != 3; ++j) {
957 const int jdx = pj * 4 + j;
958 d_value_d_dof[0](idx, jdx) +=
959 lambda_area * d_normal_d_u[i][pj][j] * weight[pi];
960 d_value_d_dof[0](dof_numbering_idx + pi * 3 + i, jdx) +=
961 -lambda_area
962 * (d_normal_d_u[i][pj][j] * c_weight[pi]
963 + normal[i] * d_c_weight[pi][pj][j]);
964 d_value_d_dof[0](
965 dof_numbering_idx + pi * 3 + i,
966 dof_numbering_idx + pj * 3 + j) +=
967 -lambda_area * normal[i] * d_c_weight[pi][pj + 3][j];
968 }
969 d_value_d_dof[0](idx, pj * 4 + 3) +=
970 area * normal[i] * weight[pj] * weight[pi];
971 d_value_d_dof[0](dof_numbering_idx + pi * 3 + i, pj * 4 + 3) +=
972 -area * normal[i] * weight[pj] * c_weight[pi];
973 }
974 }
975 for (int pj = 0; pj != 3; ++pj) {
976 for (int j = 0; j != 3; ++j) {
977 d_value_d_dof[0](pi * 4 + 3, pj * 4 + j) +=
978 -weight[pi] * d_w_d_g * d_g_d_u[pj][j];
979 d_value_d_dof[0](pi * 4 + 3, dof_numbering_idx + pj * 3 + j) +=
980 -weight[pi] * d_w_d_g * d_g_d_u[pj + 3][j];
981 }
982 d_value_d_dof[0](pi * 4 + 3, pj * 4 + 3) +=
983 -weight[pi] * d_w_d_lambda * weight[pj];
984 }
985 }
986 }
987 }
988
989 if (list_contact_element.empty()) {
990 // We add unconnected small rigidity on the contact dof.
991 for (int i = 0; i != 3; ++i) { dof_numbering[i] = dof_numbering[i * 4 + 3]; }
992 const double warea = area * composit_integration.nb_point[integration_level] / 3;
993 dof_numbering.resize(3);
994 if (!value.is_null()) {
995 value.resize(3);
996 for (int i = 0; i != 3; ++i) { value[i] = warea * dof(dof_numbering[i], 0); }
997 }
998 if (!d_value_d_dof_flags.empty() && d_value_d_dof_flags[0]) {
999 d_value_d_dof[0].resize(
1000 3, 3
1001#ifndef USE_UNSYMMETRIC
1002 ,
1003 true
1004#endif
1005 );
1006 d_value_d_dof[0].set_zero();
1007 for (int i = 0; i != 3; ++i) { d_value_d_dof[0](i, i) = warea; }
1008 }
1009 }
1010#ifndef USE_UNSYMMETRIC
1011 else if (!d_value_d_dof_flags.empty() && d_value_d_dof_flags[0])
1012 for (size_t i = 0; i != d_value_d_dof[0].size1(); ++i) {
1013 for (int j = 0; j != i; ++j) { d_value_d_dof[0](i, j) *= 0.5; }
1014 }
1015#endif
1016
1017 for (size_t i = 1; i < d_value_d_dof_flags.size(); ++i) {
1018 if (d_value_d_dof_flags[i]) {
1019 d_value_d_dof[i].resize(
1020 dof_numbering.size(), dof_numbering.size()
1021#ifndef USE_UNSYMMETRIC
1022 ,
1023 true
1024#endif
1025 );
1026 d_value_d_dof[i].set_zero();
1027 }
1028 }
1029 if (!d_value_d_time.is_null()) {
1030 d_value_d_time.resize(dof_numbering.size());
1031 d_value_d_time.set_zero();
1032 }
1033 }
1034
1035 static type_t type;
1036
1037private:
1038 // The nodes list.
1039 Node* nodes[3];
1040
1041 // The property object.
1042 ContactMechanics* property{};
1043
1044 static TriangleOnePointCompositIntegration<5> composit_integration;
1045};
1046} // namespace b2000
1047
1048#endif // B2ELEMENT_CONTACT_H_
THROW
#define THROW
Definition b2exception.H:198
b2tensor_calculus.H
b2000::Element::nonlinear
@ nonlinear
Definition b2element.H:619
b2000::Element::linear
@ linear
Definition b2element.H:615
b2000::TypedElement< double >::get_value
virtual void get_value(Model &model, const bool linear, const EquilibriumSolution equilibrium_solution, const double time, const double delta_time, const b2linalg::Matrix< double, b2linalg::Mrectangle_constref > &dof, GradientContainer *gradient_container, SolverHints *solver_hints, b2linalg::Index &dof_numbering, b2linalg::Vector< double, b2linalg::Vdense > &value, const std::vector< bool > &d_value_d_dof_flags, std::vector< b2linalg::Matrix< double, b2linalg::Mpacked > > &d_value_d_dof, b2linalg::Vector< double, b2linalg::Vdense > &d_value_d_time)
Compute the internal forces and their derivatives of the element.
Definition b2element.H:1640
b2000
Contains the base classes for implementing Finite Elements.
Definition b2boundary_condition.H:32
b2000::normalise_3
T normalise_3(T a[3])
Definition b2tensor_calculus.H:418
b2000::outer_product_3
void outer_product_3(const T1 a[3], const T2 b[3], T3 c[3])
Definition b2tensor_calculus.H:389
b2000::dot_3
T dot_3(const T a[3], const T b[3])
Definition b2tensor_calculus.H:328