v0.14.0
Loading...
Searching...
No Matches
SimpleContact.cpp
Go to the documentation of this file.
1/* \file SimpleContact.cpp
2 \brief Implementation of simple contact element
3*/
4
5/* This file is part of MoFEM.
6 * MoFEM is free software: you can redistribute it and/or modify it under
7 * the terms of the GNU Lesser General Public License as published by the
8 * Free Software Foundation, either version 3 of the License, or (at your
9 * option) any later version.
10 *
11 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14 * License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18
19using namespace MoFEM;
20using namespace boost::numeric;
21
22#include <IntegrationRules.hpp>
23
24constexpr double SimpleContactProblem::TOL;
27
31 if (newtonCotes) {
32 int rule = order + 2;
33 int nb_gauss_pts = IntRules::NCC::triangle_ncc_order_num(rule);
34 gaussPtsMaster.resize(3, nb_gauss_pts, false);
35 gaussPtsSlave.resize(3, nb_gauss_pts, false);
36 double xy_coords[2 * nb_gauss_pts];
37 double w_array[nb_gauss_pts];
38 IntRules::NCC::triangle_ncc_rule(rule, nb_gauss_pts, xy_coords, w_array);
39
40 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
41 gaussPtsMaster(0, gg) = xy_coords[gg * 2];
42 gaussPtsMaster(1, gg) = xy_coords[gg * 2 + 1];
43 gaussPtsMaster(2, gg) = w_array[gg];
44 gaussPtsSlave(0, gg) = xy_coords[gg * 2];
45 gaussPtsSlave(1, gg) = xy_coords[gg * 2 + 1];
46 gaussPtsSlave(2, gg) = w_array[gg];
47 }
48 } else {
50 }
52}
53
54template <bool CONVECT_MASTER>
58
59 auto get_material_dofs_from_coords = [&]() {
61 materialCoords.resize(18, false);
62 int num_nodes;
63 const EntityHandle *conn;
64 CHKERR fePtr->mField.get_moab().get_connectivity(fePtr->getFEEntityHandle(),
65 conn, num_nodes, true);
66 CHKERR fePtr->mField.get_moab().get_coords(conn, 6,
67 &*materialCoords.data().begin());
68 CHKERR fePtr->getNodeData(materialPositionsField, materialCoords, false);
70 };
71
72 auto get_dofs_data_for_slave_and_master = [&] {
73 slaveSpatialCoords.resize(3, 3, false);
74 slaveMaterialCoords.resize(3, 3, false);
75 masterSpatialCoords.resize(3, 3, false);
76 masterMaterialCoords.resize(3, 3, false);
77 for (size_t n = 0; n != 3; ++n) {
78 for (size_t d = 0; d != 3; ++d) {
79 masterSpatialCoords(n, d) = spatialCoords(3 * n + d);
80 slaveSpatialCoords(n, d) = spatialCoords(3 * (n + 3) + d);
81 masterMaterialCoords(n, d) = materialCoords(3 * n + d);
82 slaveMaterialCoords(n, d) = materialCoords(3 * (n + 3) + d);
83 }
84 }
85 };
86
87 auto calculate_shape_base_functions = [&](const int nb_gauss_pts) {
89 if (nb_gauss_pts != fePtr->gaussPtsMaster.size2())
90 SETERRQ2(
91 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
92 "Inconsistent size of slave and master integration points (%d != %d)",
93 nb_gauss_pts, fePtr->gaussPtsMaster.size2());
94 slaveN.resize(nb_gauss_pts, 3, false);
95 masterN.resize(nb_gauss_pts, 3, false);
96 CHKERR Tools::shapeFunMBTRI(&slaveN(0, 0), &fePtr->gaussPtsSlave(0, 0),
97 &fePtr->gaussPtsSlave(1, 0), nb_gauss_pts);
98 CHKERR Tools::shapeFunMBTRI(&masterN(0, 0), &fePtr->gaussPtsMaster(0, 0),
99 &fePtr->gaussPtsMaster(1, 0), nb_gauss_pts);
101 };
102
103 auto get_diff_ksi_master = [&]() -> MatrixDouble & {
104 if (CONVECT_MASTER)
105 return diffKsiMaster;
106 else
107 return diffKsiSlave;
108 };
109
110 auto get_diff_ksi_slave = [&]() -> MatrixDouble & {
111 if (CONVECT_MASTER)
112 return diffKsiSlave;
113 else
114 return diffKsiMaster;
115 };
116
117 auto get_slave_material_coords = [&]() -> MatrixDouble & {
118 if (CONVECT_MASTER)
119 return slaveMaterialCoords;
120 else
121 return masterMaterialCoords;
122 };
123
124 auto get_master_gauss_pts = [&]() -> MatrixDouble & {
125 if (CONVECT_MASTER)
126 return fePtr->gaussPtsMaster;
127 else
128 return fePtr->gaussPtsSlave;
129 };
130
131 auto get_slave_spatial_coords = [&]() -> MatrixDouble & {
132 if (CONVECT_MASTER)
133 return slaveSpatialCoords;
134 else
135 return masterSpatialCoords;
136 };
137
138 auto get_master_spatial_coords = [&]() -> MatrixDouble & {
139 if (CONVECT_MASTER)
140 return masterSpatialCoords;
141 else
142 return slaveSpatialCoords;
143 };
144
145 auto get_slave_n = [&]() -> MatrixDouble & {
146 if (CONVECT_MASTER)
147 return slaveN;
148 else
149 return masterN;
150 };
151
152 auto get_master_n = [&]() -> MatrixDouble & {
153 if (CONVECT_MASTER)
154 return masterN;
155 else
156 return slaveN;
157 };
158
159 auto convect_points = [get_diff_ksi_master, get_diff_ksi_slave,
160 get_slave_material_coords, get_master_gauss_pts,
161 get_slave_spatial_coords, get_master_spatial_coords,
162 get_slave_n, get_master_n](const int nb_gauss_pts) {
163 MatrixDouble3by3 A(2, 2);
164 MatrixDouble3by3 invA(2, 2);
165 VectorDouble3 F(2);
166 MatrixDouble3by3 inv_matA(2, 2);
167 VectorDouble3 copy_F(2);
169 &copy_F[1]);
171 &inv_matA(0, 0), &inv_matA(0, 1), &inv_matA(1, 0), &inv_matA(1, 1));
172
173 auto get_t_coords = [](auto &m) {
175 &m(0, 0), &m(0, 1), &m(0, 2)};
176 };
177
178 auto get_t_xi = [](auto &m) {
180 &m(1, 0)};
181 };
182
183 auto get_t_diff = []() {
186 };
187
188 auto get_t_tau = []() {
190 return t_tau;
191 };
192
193 auto get_t_x = []() {
195 return t_x;
196 };
197
198 auto get_t_F = [&]() {
200 };
201
202 auto get_t_A = [&](auto &m) {
204 &m(0, 0), &m(0, 1), &m(1, 0), &m(1, 1)};
205 };
206
207 auto get_diff_ksi = [](auto &m, const int gg) {
209 &m(0, gg), &m(1, gg), &m(2, gg), &m(3, gg), &m(4, gg), &m(5, gg));
210 };
211
217
218 get_diff_ksi_master().resize(6, 3 * nb_gauss_pts, false);
219 get_diff_ksi_slave().resize(6, 3 * nb_gauss_pts, false);
220
221 auto t_xi_master = get_t_xi(get_master_gauss_pts());
222 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
223
224 auto t_tau = get_t_tau();
225 auto t_x_slave = get_t_x();
226 auto t_x_master = get_t_x();
227 auto t_mat = get_t_A(A);
228 auto t_f = get_t_F();
229
230 auto newton_solver = [&]() {
231 auto get_values = [&]() {
232 t_tau(i, I) = 0;
233 t_x_slave(i) = 0;
234 t_x_master(i) = 0;
235
236 auto t_slave_material_coords =
237 get_t_coords(get_slave_material_coords());
238 auto t_slave_spatial_coords =
239 get_t_coords(get_slave_spatial_coords());
240 auto t_master_spatial_coords =
241 get_t_coords(get_master_spatial_coords());
242 double *slave_base = &get_slave_n()(gg, 0);
243 double *master_base = &get_master_n()(gg, 0);
244 auto t_diff = get_t_diff();
245 for (size_t n = 0; n != 3; ++n) {
246
247 t_tau(i, J) += t_diff(J) * t_slave_material_coords(i);
248 t_x_slave(i) += (*slave_base) * t_slave_spatial_coords(i);
249 t_x_master(i) += (*master_base) * t_master_spatial_coords(i);
250
251 ++t_diff;
252 ++t_slave_material_coords;
253 ++t_slave_spatial_coords;
254 ++t_master_spatial_coords;
255 ++slave_base;
256 ++master_base;
257 }
258 };
259
260 auto assemble = [&]() {
261 t_mat(I, J) = 0;
262 auto t_master_spatial_coords =
263 get_t_coords(get_master_spatial_coords());
264 auto t_diff = get_t_diff();
265 for (size_t n = 0; n != 3; ++n) {
266 t_mat(I, J) += t_diff(J) * t_tau(i, I) * t_master_spatial_coords(i);
267 ++t_diff;
268 ++t_master_spatial_coords;
269 };
270 t_f(I) = t_tau(i, I) * (t_x_slave(i) - t_x_master(i));
271 };
272
273 auto update = [&]() {
274 t_xi_master(I) += t_f(I);
275 get_master_n()(gg, 0) =
276 Tools::shapeFunMBTRI0(t_xi_master(0), t_xi_master(1));
277 get_master_n()(gg, 1) =
278 Tools::shapeFunMBTRI1(t_xi_master(0), t_xi_master(1));
279 get_master_n()(gg, 2) =
280 Tools::shapeFunMBTRI2(t_xi_master(0), t_xi_master(1));
281 };
282
283 auto invert_2_by_2 = [&](MatrixDouble3by3 &inv_mat_A,
284 MatrixDouble3by3 &mat_A) {
285 double det_A;
286 CHKERR determinantTensor2by2(mat_A, det_A);
287 CHKERR invertTensor2by2(mat_A, det_A, inv_mat_A);
288 };
289
290 auto linear_solver = [&]() {
291 invert_2_by_2(inv_matA, A);
292 t_copy_F(J) = t_f(J);
293 t_f(I) = t_inv_matA(I, J) * t_copy_F(J);
294 };
295
296 auto invert_A = [&]() { invert_2_by_2(invA, A); };
297
298 auto nonlinear_solve = [&]() {
299 constexpr double tol = 1e-12;
300 constexpr int max_it = 10;
301 int it = 0;
302 double eps;
303
304 do {
305
306 get_values();
307 assemble();
308 linear_solver();
309 update();
310
311 eps = norm_2(F);
312
313 } while (eps > tol && (it++) < max_it);
314 };
315
316 nonlinear_solve();
317 get_values();
318 assemble();
319 invert_A();
320
321 auto get_diff_slave = [&]() {
322 auto t_inv_A = get_t_A(invA);
323 auto t_diff_xi_slave = get_diff_ksi(get_diff_ksi_slave(), 3 * gg);
324 double *slave_base = &get_slave_n()(gg, 0);
325 for (size_t n = 0; n != 3; ++n) {
326 t_diff_xi_slave(I, i) = t_inv_A(I, J) * t_tau(i, J) * (*slave_base);
327 ++t_diff_xi_slave;
328 ++slave_base;
329 }
330 };
331
332 auto get_diff_master = [&]() {
333 auto t_inv_A = get_t_A(invA);
334 auto t_diff_xi_master = get_diff_ksi(get_diff_ksi_master(), 3 * gg);
335 auto t_diff = get_t_diff();
336 double *master_base = &get_master_n()(gg, 0);
338 t_diff_A(I, J, K, L) = -t_inv_A(I, K) * t_inv_A(L, J);
339 for (size_t n = 0; n != 3; ++n) {
340 t_diff_xi_master(I, i) =
341 (t_diff_A(I, J, K, L) * (t_f(J) * t_diff(L))) * t_tau(i, K) -
342 t_inv_A(I, J) * t_tau(i, J) * (*master_base);
343 ++t_diff_xi_master;
344 ++master_base;
345 ++t_diff;
346 }
347 };
348
349 get_diff_master();
350 get_diff_slave();
351 };
352
353 newton_solver();
354
355 ++t_xi_master;
356 }
357 };
358
359 const int nb_gauss_pts = fePtr->gaussPtsSlave.size2();
360 CHKERR fePtr->getNodeData(sparialPositionsField, spatialCoords);
361 CHKERR get_material_dofs_from_coords();
362 get_dofs_data_for_slave_and_master();
363 CHKERR calculate_shape_base_functions(nb_gauss_pts);
364 convect_points(nb_gauss_pts);
365
367}
368
373 CHKERR convectPtr->convectSlaveIntegrationPts<true>();
375}
376
381 CHKERR convectPtr->convectSlaveIntegrationPts<false>();
383}
385 EntityType type,
386 EntData &data) {
388
389 if (data.getFieldData().size() == 0)
391
392 if (type != MBVERTEX)
394
396
397 auto get_tensor_vec = [](VectorDouble &n) {
398 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
399 };
400
401 const double *normal_slave_ptr = &getNormalSlave()[0];
402
403 commonDataSimpleContact->normalVectorSlavePtr.get()->resize(3);
404 commonDataSimpleContact->normalVectorSlavePtr.get()->clear();
405
406 auto t_normal =
407 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
408
409 for (int ii = 0; ii != 3; ++ii)
410 t_normal(ii) = normal_slave_ptr[ii];
411
412 const double normal_length = sqrt(t_normal(i) * t_normal(i));
413 t_normal(i) = t_normal(i) / normal_length;
414
415 commonDataSimpleContact->areaSlave = 0.5 * normal_length;
416
418}
419
421 EntityType type,
422 EntData &data) {
424
425 if (data.getFieldData().size() == 0)
427
428 if (type != MBVERTEX)
430
432
433 auto get_tensor_vec = [](VectorDouble &n) {
434 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
435 };
436
437 const double *normal_master_ptr = &getNormalMaster()[0];
438 commonDataSimpleContact->normalVectorMasterPtr.get()->resize(3);
439 commonDataSimpleContact->normalVectorMasterPtr.get()->clear();
440
441 auto t_normal =
442 get_tensor_vec(commonDataSimpleContact->normalVectorMasterPtr.get()[0]);
443 for (int ii = 0; ii != 3; ++ii)
444 t_normal(ii) = normal_master_ptr[ii];
445
446 const double normal_length = sqrt(t_normal(i) * t_normal(i));
447 t_normal(i) = t_normal(i) / normal_length;
448 commonDataSimpleContact->areaMaster = 0.5 * normal_length;
449
451}
452
454 int side, EntityType type, EntData &data) {
456 const int nb_dofs = data.getFieldData().size();
457
458 if (nb_dofs == 0)
460
461 const int nb_gauss_pts = data.getN().size1();
462
463 if (type == MBVERTEX) {
464 commonDataSimpleContact->positionAtGaussPtsMasterPtr.get()->resize(
465 3, nb_gauss_pts, false);
466
467 commonDataSimpleContact->positionAtGaussPtsMasterPtr.get()->clear();
468 }
469
470 auto t_position_master = getFTensor1FromMat<3>(
471 *commonDataSimpleContact->positionAtGaussPtsMasterPtr);
472
473 int nb_base_fun_col = data.getFieldData().size() / 3;
474
476
477 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
478 FTensor::Tensor0<double *> t_base_master(&data.getN()(gg, 0));
479
480 FTensor::Tensor1<double *, 3> t_field_data_master(
481 &data.getFieldData()[0], &data.getFieldData()[1],
482 &data.getFieldData()[2], 3);
483
484 for (int bb = 0; bb != nb_base_fun_col; ++bb) {
485 t_position_master(i) += t_base_master * t_field_data_master(i);
486
487 ++t_base_master;
488 ++t_field_data_master;
489 }
490 ++t_position_master;
491 }
492
494}
495
498 int side, EntityType type, EntData &data) {
500 const int nb_dofs = data.getFieldData().size();
501
502 if (nb_dofs == 0)
504
505 const int nb_gauss_pts = data.getN().size1();
506
507 auto t_new_spat_pos_master = getFTensor1FromMat<3>(
508 *commonDataSimpleContact->positionAtGaussPtsMasterPtr);
509
510 int nb_base_fun_col = data.getFieldData().size() / 3;
511
513
514 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
515 FTensor::Tensor0<double *> t_base_master(&data.getN()(gg, 0));
516
517 FTensor::Tensor1<double *, 3> t_field_data_master(
518 &data.getFieldData()[0], &data.getFieldData()[1],
519 &data.getFieldData()[2], 3);
520
521 for (int bb = 0; bb != nb_base_fun_col; ++bb) {
522 t_new_spat_pos_master(i) -= t_base_master * t_field_data_master(i);
523
524 ++t_base_master;
525 ++t_field_data_master;
526 }
527 ++t_new_spat_pos_master;
528 }
529
531}
532
535 int side, EntityType type, EntData &data) {
537 const int nb_dofs = data.getFieldData().size();
538
539 if (nb_dofs == 0)
541
542 const int nb_gauss_pts = data.getN().size1();
543
544 auto t_new_spat_pos_master = getFTensor1FromMat<3>(
545 *commonDataSimpleContact->positionAtGaussPtsMasterPtr);
546
547 int nb_base_fun_col = data.getFieldData().size() / 3;
548
550
551 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
552 FTensor::Tensor0<double *> t_base_master(&data.getN()(gg, 0));
553
554 FTensor::Tensor1<double *, 3> t_field_data_master(
555 &data.getFieldData()[0], &data.getFieldData()[1],
556 &data.getFieldData()[2], 3);
557
558 for (int bb = 0; bb != nb_base_fun_col; ++bb) {
559 t_new_spat_pos_master(i) += t_base_master * t_field_data_master(i);
560
561 ++t_base_master;
562 ++t_field_data_master;
563 }
564 ++t_new_spat_pos_master;
565 }
566
568}
569
571 int side, EntityType type, EntData &data) {
573 const int nb_dofs = data.getFieldData().size();
574
575 if (nb_dofs == 0)
577
578 const int nb_gauss_pts = data.getN().size1();
579
580 if (type == MBVERTEX) {
581 commonDataSimpleContact->positionAtGaussPtsSlavePtr.get()->resize(
582 3, nb_gauss_pts, false);
583
584 commonDataSimpleContact->positionAtGaussPtsSlavePtr.get()->clear();
585 }
586
587 auto t_position_slave = getFTensor1FromMat<3>(
588 *commonDataSimpleContact->positionAtGaussPtsSlavePtr);
589
590 int nb_base_fun_col = data.getFieldData().size() / 3;
591
593
594 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
595 FTensor::Tensor0<double *> t_base_slave(&data.getN()(gg, 0));
596
597 FTensor::Tensor1<double *, 3> t_field_data_slave(
598 &data.getFieldData()[0], &data.getFieldData()[1],
599 &data.getFieldData()[2], 3);
600
601 for (int bb = 0; bb != nb_base_fun_col; ++bb) {
602 t_position_slave(i) += t_base_slave * t_field_data_slave(i);
603
604 ++t_base_slave;
605 ++t_field_data_slave;
606 }
607 ++t_position_slave;
608 }
609
611}
612
614 int side, EntityType type, EntData &data) {
616 const int nb_dofs = data.getFieldData().size();
617
618 if (nb_dofs == 0)
620
621 const int nb_gauss_pts = data.getN().size1();
622
623 auto t_new_spat_pos_slave = getFTensor1FromMat<3>(
624 *commonDataSimpleContact->positionAtGaussPtsSlavePtr);
625
626 int nb_base_fun_col = data.getFieldData().size() / 3;
627
629
630 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
631 FTensor::Tensor0<double *> t_base_slave(&data.getN()(gg, 0));
632
633 FTensor::Tensor1<double *, 3> t_field_data_slave(
634 &data.getFieldData()[0], &data.getFieldData()[1],
635 &data.getFieldData()[2], 3);
636
637 for (int bb = 0; bb != nb_base_fun_col; ++bb) {
638 t_new_spat_pos_slave(i) -= t_base_slave * t_field_data_slave(i);
639
640 ++t_base_slave;
641 ++t_field_data_slave;
642 }
643 ++t_new_spat_pos_slave;
644 }
645
647}
648
651 int side, EntityType type, EntData &data) {
653 const int nb_dofs = data.getFieldData().size();
654
655 if (nb_dofs == 0)
657
658 const int nb_gauss_pts = data.getN().size1();
659
660 auto t_new_spat_pos_slave = getFTensor1FromMat<3>(
661 *commonDataSimpleContact->positionAtGaussPtsSlavePtr);
662
663 int nb_base_fun_col = data.getFieldData().size() / 3;
664
666
667 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
668 FTensor::Tensor0<double *> t_base_slave(&data.getN()(gg, 0));
669
670 FTensor::Tensor1<double *, 3> t_field_data_slave(
671 &data.getFieldData()[0], &data.getFieldData()[1],
672 &data.getFieldData()[2], 3);
673
674 for (int bb = 0; bb != nb_base_fun_col; ++bb) {
675 t_new_spat_pos_slave(i) += t_base_slave * t_field_data_slave(i);
676
677 ++t_base_slave;
678 ++t_field_data_slave;
679 }
680 ++t_new_spat_pos_slave;
681 }
682
684}
685
687 EntityType type,
688 EntData &data) {
690
691 if (data.getFieldData().size() == 0)
693
694 if (type != MBVERTEX)
696
697 auto get_tensor_vec = [](VectorDouble &n) {
698 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
699 };
700
701 const int nb_gauss_pts = data.getN().size1();
702
703 commonDataSimpleContact->gapPtr.get()->resize(nb_gauss_pts);
704 commonDataSimpleContact->gapPtr.get()->clear();
705
707
708 auto t_position_master_gp = getFTensor1FromMat<3>(
709 *commonDataSimpleContact->positionAtGaussPtsMasterPtr);
710 auto t_position_slave_gp = getFTensor1FromMat<3>(
711 *commonDataSimpleContact->positionAtGaussPtsSlavePtr);
712
713 auto t_gap_ptr = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
714
715 auto t_normal_at_gp =
716 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
717
718 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
719 t_gap_ptr -=
720 t_normal_at_gp(i) * (t_position_slave_gp(i) - t_position_master_gp(i));
721 ++t_position_slave_gp;
722 ++t_position_master_gp;
723 ++t_gap_ptr;
724 } // for gauss points
725
727}
728
730 int side, EntityType type, EntData &data) {
732
733 if (data.getFieldData().size() == 0)
735
736 const int nb_gauss_pts = data.getN().size1();
737
738 if (type == MBVERTEX) {
739 commonDataSimpleContact->lagMultAtGaussPtsPtr.get()->resize(nb_gauss_pts);
740 commonDataSimpleContact->lagMultAtGaussPtsPtr.get()->clear();
741 }
742
743 int nb_base_fun_row = data.getFieldData().size();
744
745 auto t_lagrange_slave =
746 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
747
748 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
749 FTensor::Tensor0<double *> t_base_lambda(&data.getN()(gg, 0));
750
751 FTensor::Tensor0<double *> t_field_data_slave(&data.getFieldData()[0]);
752 for (int bb = 0; bb != nb_base_fun_row; ++bb) {
753 t_lagrange_slave += t_base_lambda * t_field_data_slave;
754 ++t_base_lambda;
755 ++t_field_data_slave;
756 }
757 ++t_lagrange_slave;
758 }
759
761}
762
764 int side, EntityType type, EntData &data) {
766
767 if (data.getFieldData().size() == 0)
769
770 if (type != MBVERTEX)
772
773 auto get_tensor_vec = [](VectorDouble &n) {
774 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
775 };
776
777 const int nb_gauss_pts = data.getN().size1();
778
779 commonDataSimpleContact->augmentedLambdasPtr.get()->resize(nb_gauss_pts);
780 commonDataSimpleContact->augmentedLambdasPtr.get()->clear();
781
783
784 auto t_aug_lambda_ptr =
785 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
786
787 auto t_gap_ptr = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
788
789 auto t_lagrange_slave =
790 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
791
792 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
793 t_aug_lambda_ptr += t_lagrange_slave + cN * t_gap_ptr;
794 ++t_aug_lambda_ptr;
795 ++t_lagrange_slave;
796 ++t_gap_ptr;
797 } // for gauss points
798
800}
801
804 doWork(int row_side, int col_side, EntityType row_type, EntityType col_type,
805 EntData &row_data, EntData &col_data) {
807
808 // Both sides are needed since both sides contribute their shape
809 // function to the stiffness matrix
810 const int nb_row = row_data.getIndices().size();
811 const int nb_col = col_data.getIndices().size();
812
813 if (nb_row && nb_col) {
814 const int nb_gauss_pts = row_data.getN().size1();
815
816 int nb_base_fun_row = row_data.getFieldData().size() / 3;
817 int nb_base_fun_col = col_data.getFieldData().size();
818
819 const double area_master =
820 commonDataSimpleContact->areaSlave; // same area in master and slave
821
822 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
823 return FTensor::Tensor1<double *, 3>(&m(r + 0, c), &m(r + 1, c),
824 &m(r + 2, c));
825 };
826
827 auto get_tensor_vec = [](VectorDouble &n) {
828 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
829 };
830
832
833 NN.resize(3 * nb_base_fun_row, nb_base_fun_col, false);
834 NN.clear();
835
836 auto const_unit_n =
837 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr.get()));
838
839 auto t_w = getFTensor0IntegrationWeightSlave();
840
841 auto t_aug_lambda_ptr =
842 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
843
844 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
845
846 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
847 ++t_w;
848 ++t_aug_lambda_ptr;
849 continue;
850 }
851 double val_m = t_w * area_master;
852 auto t_base_master = row_data.getFTensor0N(gg, 0);
853
854 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
855 auto t_assemble_m = get_tensor_from_mat(NN, 3 * bbr, 0);
856 auto t_base_lambda = col_data.getFTensor0N(gg, 0);
857 const double m = val_m * t_base_master;
858 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
859 const double n = m * t_base_lambda;
860 t_assemble_m(i) += n * const_unit_n(i);
861 ++t_assemble_m;
862 ++t_base_lambda; // update cols slave
863 }
864 ++t_base_master; // update rows master
865 }
866 ++t_w;
867 ++t_aug_lambda_ptr;
868 }
869
870 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
871 ADD_VALUES);
872 }
873
875}
876
879 int row_side, int col_side, EntityType row_type, EntityType col_type,
880 EntData &row_data, EntData &col_data) {
882
883 // Both sides are needed since both sides contribute their shape
884 // function to the stiffness matrix
885 const int nb_row = row_data.getIndices().size();
886 const int nb_col = col_data.getIndices().size();
887
888 if (nb_row && nb_col) {
889 const int nb_gauss_pts = row_data.getN().size1();
890
891 int nb_base_fun_row = row_data.getFieldData().size() / 3;
892 int nb_base_fun_col = col_data.getFieldData().size();
893
894 const double area_slave =
895 commonDataSimpleContact->areaSlave; // same area in master and slave
896
897 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
898 return FTensor::Tensor1<double *, 3>(&m(r + 0, c), &m(r + 1, c),
899 &m(r + 2, c));
900 };
901
902 auto get_tensor_vec = [](VectorDouble &n) {
903 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
904 };
905
907
908 NN.resize(3 * nb_base_fun_row, nb_base_fun_col, false);
909 NN.clear();
910
911 auto const_unit_n =
912 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr.get()));
913
914 auto t_w = getFTensor0IntegrationWeightSlave();
915
916 auto t_aug_lambda_ptr =
917 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
918
919 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
920
921 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
922 ++t_w;
923 ++t_aug_lambda_ptr;
924 continue;
925 }
926
927 double val_m = t_w * area_slave;
928 auto t_base_master = row_data.getFTensor0N(gg, 0);
929
930 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
931 auto t_assemble_m = get_tensor_from_mat(NN, 3 * bbr, 0);
932 auto t_base_lambda = col_data.getFTensor0N(gg, 0);
933 const double m = val_m * t_base_master;
934 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
935 const double n = m * t_base_lambda;
936 t_assemble_m(i) -= n * const_unit_n(i);
937
938 ++t_assemble_m;
939 ++t_base_lambda; // update cols slave
940 }
941 ++t_base_master; // update rows master
942 }
943 ++t_w;
944 ++t_aug_lambda_ptr;
945 }
946
947 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
948 ADD_VALUES);
949 }
951}
952
954 int row_side, int col_side, EntityType row_type, EntityType col_type,
955 EntData &row_data, EntData &col_data) {
957
958 const int nb_row = row_data.getIndices().size();
959 const int nb_col = col_data.getIndices().size();
960
961 if (nb_row && nb_col) {
962 const int nb_gauss_pts = row_data.getN().size1();
963 const double area_slave =
964 commonDataSimpleContact->areaSlave; // same area in master and slave
965
966 NN.resize(nb_row, nb_col, false);
967 NN.clear();
968
969 auto t_lagrange_slave =
970 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
971 auto t_w = getFTensor0IntegrationWeightSlave();
972
973 auto t_aug_lambda_ptr =
974 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
975
976 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
977
978 if (t_aug_lambda_ptr <= 0) {
979 ++t_w;
980 ++t_aug_lambda_ptr;
981 ++t_lagrange_slave;
982 continue;
983 }
984
985 const double val_s = -t_w * area_slave / cN;
986
988 &*NN.data().begin());
989
990 auto t_base_lambda_row = row_data.getFTensor0N(gg, 0);
991 for (int bbr = 0; bbr != nb_row; ++bbr) {
992 auto t_base_lambda_col = col_data.getFTensor0N(gg, 0);
993 const double s = val_s * t_base_lambda_row;
994 for (int bbc = 0; bbc != nb_col; ++bbc) {
995
996 t_mat += s * t_base_lambda_col;
997
998 ++t_mat;
999 ++t_base_lambda_col; // update cols
1000 }
1001 ++t_base_lambda_row; // update rows
1002 }
1003
1004 ++t_lagrange_slave;
1005 ++t_w;
1006 ++t_aug_lambda_ptr;
1007 }
1008 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1009 ADD_VALUES);
1010 }
1011
1013}
1014
1017 int row_side, int col_side, EntityType row_type, EntityType col_type,
1018 EntData &row_data, EntData &col_data) {
1020
1021 const int nb_row = row_data.getIndices().size();
1022 const int nb_col = col_data.getIndices().size();
1023
1024 if (nb_row && nb_col) {
1025
1026 const int nb_gauss_pts = row_data.getN().size1();
1027 int nb_base_fun_row = row_data.getFieldData().size();
1028 int nb_base_fun_col = col_data.getFieldData().size() / 3;
1029
1030 const double area_master =
1031 commonDataSimpleContact->areaSlave; // same area in master and slave
1032
1033 const double area_slave =
1034 commonDataSimpleContact->areaSlave; // same area in master and slave
1035
1036 NN.resize(nb_base_fun_row, 3 * nb_base_fun_col, false);
1037 NN.clear();
1038
1039 auto get_tensor_from_vec = [](VectorDouble &n) {
1040 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
1041 };
1042
1044
1045 auto t_const_unit_n =
1046 get_tensor_from_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
1047
1048 auto t_const_unit_master =
1049 get_tensor_from_vec(*(commonDataSimpleContact->normalVectorMasterPtr));
1050
1051 auto t_aug_lambda_ptr =
1052 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1053
1054 auto t_w = getFTensor0IntegrationWeightSlave();
1055 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1056
1057 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1058 ++t_w;
1059 ++t_aug_lambda_ptr;
1060 continue;
1061 }
1062
1063 const double val_m = t_w * area_slave;
1064
1065 auto t_base_lambda = row_data.getFTensor0N(gg, 0);
1066 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
1067
1068 auto t_base_master = col_data.getFTensor0N(gg, 0);
1070 &NN(bbr, 0), &NN(bbr, 1), &NN(bbr, 2)};
1071 const double m = val_m * t_base_lambda;
1072
1073 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
1074
1075 t_mat(i) += t_const_unit_n(i) * m * t_base_master;
1076
1077 ++t_base_master; // update rows
1078 ++t_mat;
1079 }
1080 ++t_base_lambda; // update cols master
1081 }
1082
1083 ++t_aug_lambda_ptr;
1084 ++t_w;
1085 }
1086
1087 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1088 ADD_VALUES);
1089 }
1090
1092}
1093
1096 int row_side, int col_side, EntityType row_type, EntityType col_type,
1097 EntData &row_data, EntData &col_data) {
1099
1100 const int nb_row = row_data.getIndices().size();
1101 const int nb_col = col_data.getIndices().size();
1102
1103 if (nb_row && nb_col) {
1104
1105 const int nb_gauss_pts = row_data.getN().size1();
1106 int nb_base_fun_row = row_data.getFieldData().size();
1107 int nb_base_fun_col = col_data.getFieldData().size() / 3;
1108
1109 const double area_slave =
1110 commonDataSimpleContact->areaSlave; // same area in master and slave
1111
1112 NN.resize(nb_base_fun_row, 3 * nb_base_fun_col, false);
1113 NN.clear();
1114
1115 auto get_tensor_from_vec = [](VectorDouble &n) {
1116 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
1117 };
1118
1120
1121 auto t_const_unit_n =
1122 get_tensor_from_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
1123
1124 auto t_aug_lambda_ptr =
1125 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1126
1127 auto t_w = getFTensor0IntegrationWeightSlave();
1128 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1129
1130 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1131 ++t_w;
1132 ++t_aug_lambda_ptr;
1133 continue;
1134 }
1135
1136 const double val_m = t_w * area_slave;
1137
1138 auto t_base_lambda = row_data.getFTensor0N(gg, 0);
1139 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
1140
1141 auto t_base_slave = col_data.getFTensor0N(gg, 0);
1143 &NN(bbr, 0), &NN(bbr, 1), &NN(bbr, 2)};
1144 const double m = val_m * t_base_lambda;
1145
1146 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
1147
1148 t_mat(i) -= t_const_unit_n(i) * m * t_base_slave;
1149
1150 ++t_base_slave; // update rows
1151 ++t_mat;
1152 }
1153 ++t_base_lambda; // update cols master
1154 }
1155
1156 ++t_aug_lambda_ptr;
1157 ++t_w;
1158 }
1159
1160 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1161 ADD_VALUES);
1162 }
1163
1165}
1166
1169 doWork(int row_side, int col_side, EntityType row_type, EntityType col_type,
1171 EntitiesFieldData::EntData &col_data) {
1173
1174 // Both sides are needed since both sides contribute their shape
1175 // function to the stiffness matrix
1176 const int nb_row = row_data.getIndices().size();
1177 if (!nb_row)
1179 const int nb_col = col_data.getIndices().size();
1180 if (!nb_col)
1182 const int nb_gauss_pts = row_data.getN().size1();
1183
1184 int nb_base_fun_row = row_data.getFieldData().size() / 3;
1185 int nb_base_fun_col = col_data.getFieldData().size() / 3;
1186
1187 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
1189 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1190 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1191 &m(r + 2, c + 2));
1192 };
1193
1194 auto get_tensor_vec = [](VectorDouble &n) {
1195 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
1196 };
1197
1200
1201 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
1202 NN.clear();
1203
1204 auto t_w = getFTensor0IntegrationWeightSlave();
1205
1206 const double area_master = commonDataSimpleContact->areaSlave;
1207
1208 auto normal =
1209 get_tensor_vec(commonDataSimpleContact->normalVectorSlavePtr.get()[0]);
1210
1211 auto t_aug_lambda_ptr =
1212 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1213
1214 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1215
1216 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1217 ++t_w;
1218 ++t_aug_lambda_ptr;
1219 continue;
1220 }
1221
1222 const double val_m = t_w * area_master * cN;
1223
1224 FTensor::Tensor0<double *> t_base_master_col(&col_data.getN()(gg, 0));
1225
1226 for (int bbc = 0; bbc != nb_base_fun_col; bbc++) {
1227
1228 FTensor::Tensor0<double *> t_base_master_row(&row_data.getN()(gg, 0));
1229 const double m = val_m * t_base_master_col;
1230
1231 for (int bbr = 0; bbr != nb_base_fun_row; bbr++) {
1232
1233 auto t_assemble_s = get_tensor_from_mat(NN, 3 * bbr, 3 * bbc);
1234
1235 t_assemble_s(i, j) += m * normal(i) * normal(j) * t_base_master_row;
1236
1237 ++t_base_master_row; // update rows
1238 }
1239 ++t_base_master_col; // update cols slave
1240 }
1241 ++t_w;
1242 ++t_aug_lambda_ptr;
1243 }
1244
1245 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1246 ADD_VALUES);
1247
1249}
1250
1253 doWork(int row_side, int col_side, EntityType row_type, EntityType col_type,
1255 EntitiesFieldData::EntData &col_data) {
1257
1258 // Both sides are needed since both sides contribute their shape
1259 // function to the stiffness matrix
1260 const int nb_row = row_data.getIndices().size();
1261 if (!nb_row)
1263 const int nb_col = col_data.getIndices().size();
1264 if (!nb_col)
1266 const int nb_gauss_pts = row_data.getN().size1();
1267
1268 int nb_base_fun_row = row_data.getFieldData().size() / 3;
1269 int nb_base_fun_col = col_data.getFieldData().size() / 3;
1270
1271 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
1273 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1274 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1275 &m(r + 2, c + 2));
1276 };
1277
1278 auto get_tensor_vec = [](VectorDouble &n) {
1279 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
1280 };
1281
1284
1285 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
1286 NN.clear();
1287
1288 auto t_w = getFTensor0IntegrationWeightSlave();
1289
1290 const double area_master = commonDataSimpleContact->areaSlave;
1291
1292 auto normal =
1293 get_tensor_vec(commonDataSimpleContact->normalVectorSlavePtr.get()[0]);
1294
1295 auto t_aug_lambda_ptr =
1296 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1297
1298 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1299
1300 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1301 ++t_w;
1302 ++t_aug_lambda_ptr;
1303 continue;
1304 }
1305
1306 const double val_m = t_w * area_master * cN;
1307
1308 FTensor::Tensor0<double *> t_base_slave_col(&col_data.getN()(gg, 0));
1309
1310 for (int bbc = 0; bbc != nb_base_fun_col; bbc++) {
1311
1312 FTensor::Tensor0<double *> t_base_master_row(&row_data.getN()(gg, 0));
1313 const double m = val_m * t_base_slave_col;
1314
1315 for (int bbr = 0; bbr != nb_base_fun_row; bbr++) {
1316
1317 auto t_assemble_s = get_tensor_from_mat(NN, 3 * bbr, 3 * bbc);
1318
1319 t_assemble_s(i, j) -= m * normal(i) * normal(j) * t_base_master_row;
1320
1321 ++t_base_master_row; // update rows
1322 }
1323 ++t_base_slave_col; // update cols slave
1324 }
1325 ++t_w;
1326 ++t_aug_lambda_ptr;
1327 }
1328
1329 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1330 ADD_VALUES);
1331
1333}
1334
1337 doWork(int row_side, int col_side, EntityType row_type, EntityType col_type,
1339 EntitiesFieldData::EntData &col_data) {
1341
1342 // Both sides are needed since both sides contribute their shape
1343 // function to the stiffness matrix
1344 const int nb_row = row_data.getIndices().size();
1345 if (!nb_row)
1347 const int nb_col = col_data.getIndices().size();
1348 if (!nb_col)
1350 const int nb_gauss_pts = row_data.getN().size1();
1351
1352 int nb_base_fun_row = row_data.getFieldData().size() / 3;
1353 int nb_base_fun_col = col_data.getFieldData().size() / 3;
1354
1355 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
1357 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1358 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1359 &m(r + 2, c + 2));
1360 };
1361
1362 auto get_tensor_vec = [](VectorDouble &n) {
1363 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
1364 };
1365
1368
1369 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
1370 NN.clear();
1371
1372 auto t_w = getFTensor0IntegrationWeightSlave();
1373
1374 const double area_slave = commonDataSimpleContact->areaSlave;
1375
1376 auto normal =
1377 get_tensor_vec(commonDataSimpleContact->normalVectorSlavePtr.get()[0]);
1378
1379 auto t_aug_lambda_ptr =
1380 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1381
1382 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1383
1384 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1385 ++t_w;
1386 ++t_aug_lambda_ptr;
1387 continue;
1388 }
1389
1390 const double val_s = t_w * area_slave * cN;
1391
1392 FTensor::Tensor0<double *> t_base_slave_col(&col_data.getN()(gg, 0));
1393
1394 for (int bbc = 0; bbc != nb_base_fun_col; bbc++) {
1395
1396 FTensor::Tensor0<double *> t_base_slave_row(&row_data.getN()(gg, 0));
1397 const double s = val_s * t_base_slave_col;
1398
1399 for (int bbr = 0; bbr != nb_base_fun_row; bbr++) {
1400
1401 auto t_assemble_s = get_tensor_from_mat(NN, 3 * bbr, 3 * bbc);
1402
1403 t_assemble_s(i, j) += s * normal(i) * normal(j) * t_base_slave_row;
1404
1405 ++t_base_slave_row; // update rows
1406 }
1407 ++t_base_slave_col; // update cols slave
1408 }
1409 ++t_w;
1410 ++t_aug_lambda_ptr;
1411 }
1412
1413 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1414 ADD_VALUES);
1415
1417}
1418
1421 doWork(int row_side, int col_side, EntityType row_type, EntityType col_type,
1423 EntitiesFieldData::EntData &col_data) {
1425
1426 // Both sides are needed since both sides contribute their shape
1427 // function to the stiffness matrix
1428 const int nb_row = row_data.getIndices().size();
1429 if (!nb_row)
1431 const int nb_col = col_data.getIndices().size();
1432 if (!nb_col)
1434 const int nb_gauss_pts = row_data.getN().size1();
1435
1436 int nb_base_fun_row = row_data.getFieldData().size() / 3;
1437 int nb_base_fun_col = col_data.getFieldData().size() / 3;
1438
1439 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
1441 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1442 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1443 &m(r + 2, c + 2));
1444 };
1445
1446 auto get_tensor_vec = [](VectorDouble &n) {
1447 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
1448 };
1449
1452
1453 NN.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
1454 NN.clear();
1455
1456 auto t_w = getFTensor0IntegrationWeightSlave();
1457
1458 const double area_slave = commonDataSimpleContact->areaSlave;
1459
1460 auto normal =
1461 get_tensor_vec(commonDataSimpleContact->normalVectorSlavePtr.get()[0]);
1462
1463 auto t_aug_lambda_ptr =
1464 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1465
1466 for (int gg = 0; gg != nb_gauss_pts; gg++) {
1467
1468 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1469 ++t_w;
1470 ++t_aug_lambda_ptr;
1471 continue;
1472 }
1473
1474 const double val_s = t_w * area_slave * cN;
1475
1476 FTensor::Tensor0<double *> t_base_master_col(&col_data.getN()(gg, 0));
1477
1478 for (int bbc = 0; bbc != nb_base_fun_col; bbc++) {
1479
1480 FTensor::Tensor0<double *> t_base_slave_row(&row_data.getN()(gg, 0));
1481 const double s = val_s * t_base_master_col;
1482
1483 for (int bbr = 0; bbr != nb_base_fun_row; bbr++) {
1484
1485 auto t_assemble_s = get_tensor_from_mat(NN, 3 * bbr, 3 * bbc);
1486
1487 t_assemble_s(i, j) -= s * normal(i) * normal(j) * t_base_slave_row;
1488
1489 ++t_base_slave_row; // update rows
1490 }
1491 ++t_base_master_col; // update cols slave
1492 }
1493 ++t_w;
1494 ++t_aug_lambda_ptr;
1495 }
1496
1497 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1498 ADD_VALUES);
1499
1501}
1502
1504 int side, EntityType type, EntData &data) {
1506
1507 if (type != MBVERTEX)
1509
1510 const int nb_gauss_pts = data.getN().size1();
1511
1512 commonDataSimpleContact->lagGapProdPtr.get()->resize(nb_gauss_pts);
1513 commonDataSimpleContact->lagGapProdPtr.get()->clear();
1514
1515 int nb_base_fun_row = data.getFieldData().size();
1516
1517 auto t_lagrange_slave =
1518 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
1519
1520 auto t_lag_gap_prod_slave =
1521 getFTensor0FromVec(*commonDataSimpleContact->lagGapProdPtr);
1522
1523 auto t_gap_ptr = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
1524
1525 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1526 t_lag_gap_prod_slave += t_gap_ptr * t_lagrange_slave;
1527 ++t_gap_ptr;
1528 ++t_lag_gap_prod_slave;
1529 ++t_lagrange_slave;
1530 }
1531
1533}
1534
1536 int side, EntityType type, EntData &data) {
1538
1539 const int nb_dofs = data.getIndices().size();
1540 if (nb_dofs) {
1541
1542 const int nb_gauss_pts = data.getN().size1();
1543 int nb_base_fun_col = nb_dofs / 3;
1544
1545 vecF.resize(nb_dofs, false);
1546 vecF.clear();
1547
1548 const double area_m =
1549 commonDataSimpleContact->areaMaster; // same area in master and slave
1550
1551 auto get_tensor_vec = [](VectorDouble &n, const int r) {
1552 return FTensor::Tensor1<double *, 3>(&n(r + 0), &n(r + 1), &n(r + 2));
1553 };
1554
1556
1557 auto t_lagrange_slave =
1558 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
1559
1560 auto t_const_unit_n = get_tensor_vec(
1561 commonDataSimpleContact->normalVectorSlavePtr.get()[0], 0);
1562
1563 auto t_w = getFTensor0IntegrationWeightMaster();
1564
1565 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1566
1567 double val_m = t_w * area_m;
1568
1569 auto t_base_master = data.getFTensor0N(gg, 0);
1571 &vecF[0], &vecF[1], &vecF[2]};
1572
1573 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
1574 const double m = val_m * t_base_master * t_lagrange_slave;
1575 t_assemble_m(i) -= m * t_const_unit_n(i);
1576 ++t_base_master;
1577 ++t_assemble_m;
1578 }
1579
1580 ++t_lagrange_slave;
1581 ++t_w;
1582 } // for gauss points
1583
1584 CHKERR VecSetValues(getSNESf(), data, &*vecF.begin(), ADD_VALUES);
1585 }
1587}
1588
1590 int side, EntityType type, EntData &data) {
1592
1593 const int nb_dofs = data.getIndices().size();
1594 if (nb_dofs) {
1595
1596 const int nb_gauss_pts = data.getN().size1();
1597 int nb_base_fun_col = nb_dofs / 3;
1598
1599 vecF.resize(nb_dofs, false);
1600 vecF.clear();
1601
1602 const double area_m = commonDataSimpleContact->areaSlave;
1603
1604 auto get_tensor_vec = [](VectorDouble &n, const int r) {
1605 return FTensor::Tensor1<double *, 3>(&n(r + 0), &n(r + 1), &n(r + 2));
1606 };
1607
1609
1610 auto t_lagrange_slave =
1611 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
1612
1613 auto t_const_unit_n = get_tensor_vec(
1614 commonDataSimpleContact->normalVectorSlavePtr.get()[0], 0);
1615
1616 auto t_w = getFTensor0IntegrationWeightSlave();
1617
1618 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1619
1620 double val_m = t_w * area_m;
1621
1622 auto t_base_slave = data.getFTensor0N(gg, 0);
1624 &vecF[0], &vecF[1], &vecF[2]};
1625
1626 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
1627 const double m = val_m * t_base_slave * t_lagrange_slave;
1628 t_assemble_m(i) += m * t_const_unit_n(i);
1629 ++t_base_slave;
1630 ++t_assemble_m;
1631 }
1632
1633 ++t_lagrange_slave;
1634 ++t_w;
1635 } // for gauss points
1636
1637 CHKERR VecSetValues(getSNESf(), data, &*vecF.begin(), ADD_VALUES);
1638 }
1640}
1641
1644 EntData &data) {
1646
1647 if (data.getIndices().size() == 0)
1649
1650 const int nb_gauss_pts = data.getN().size1();
1651
1652 int nb_base_fun_col = data.getFieldData().size();
1653 const double area_s =
1654 commonDataSimpleContact->areaSlave; // same area in master and slave
1655
1656 vecR.resize(nb_base_fun_col, false);
1657 vecR.clear();
1658
1659 auto t_lagrange_slave =
1660 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
1661 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
1662 auto t_w = getFTensor0IntegrationWeightSlave();
1663 const double cn_value = *cNPtr.get();
1664 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1665 const double val_s = t_w * area_s *
1667 cn_value, t_gap_gp, t_lagrange_slave);
1668 auto t_base_lambda = data.getFTensor0N(gg, 0);
1669 for (int bbr = 0; bbr != nb_base_fun_col; ++bbr) {
1670 vecR[bbr] += val_s * t_base_lambda;
1671
1672 ++t_base_lambda; // update rows
1673 }
1674
1675 ++t_lagrange_slave;
1676 ++t_gap_gp;
1677 ++t_w;
1678 } // for gauss points
1679
1680 CHKERR VecSetValues(getSNESf(), data, &vecR[0], ADD_VALUES);
1681
1683}
1684
1686 int side, EntityType type, EntData &data) {
1688
1689 if (data.getIndices().size() == 0)
1691
1692 const int nb_gauss_pts = data.getN().size1();
1693
1694 int nb_base_fun_col = data.getFieldData().size();
1695 const double area_s =
1696 commonDataSimpleContact->areaSlave; // same area in master and slave
1697
1698 vecR.resize(nb_base_fun_col, false);
1699 vecR.clear();
1700
1701 auto t_lagrange_slave =
1702 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
1703 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
1704 auto t_w = getFTensor0IntegrationWeightSlave();
1705
1706 auto t_aug_lambda_ptr =
1707 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1708
1709 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1710
1711 double branch_gg;
1712 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1713 branch_gg = -t_lagrange_slave / cN;
1714 } else {
1715 branch_gg = t_gap_gp;
1716 }
1717
1718 const double val_s = t_w * area_s * branch_gg;
1719 auto t_base_lambda = data.getFTensor0N(gg, 0);
1720 for (int bbr = 0; bbr != nb_base_fun_col; ++bbr) {
1721 vecR[bbr] += val_s * t_base_lambda;
1722
1723 ++t_base_lambda; // update rows
1724 }
1725
1726 ++t_lagrange_slave;
1727 ++t_gap_gp;
1728 ++t_aug_lambda_ptr;
1729 ++t_w;
1730 } // for gauss points
1731
1732 CHKERR VecSetValues(getSNESf(), data, &vecR[0], ADD_VALUES);
1733
1735}
1736
1739 int row_side, int col_side, EntityType row_type, EntityType col_type,
1740 EntData &row_data, EntData &col_data) {
1742
1743 // Both sides are needed since both sides contribute their shape
1744 // function to the stiffness matrix
1745 const int nb_row = row_data.getIndices().size();
1746 const int nb_col = col_data.getIndices().size();
1747
1748 if (nb_row && nb_col) {
1749 const int nb_gauss_pts = row_data.getN().size1();
1750
1751 int nb_base_fun_row = row_data.getFieldData().size() / 3;
1752 int nb_base_fun_col = col_data.getFieldData().size();
1753
1754 const double area_master =
1755 commonDataSimpleContact->areaMaster; // same area in master and slave
1756
1757 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
1758 return FTensor::Tensor1<double *, 3>(&m(r + 0, c), &m(r + 1, c),
1759 &m(r + 2, c));
1760 };
1761
1762 auto get_tensor_vec = [](VectorDouble &n) {
1763 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
1764 };
1765
1767
1768 NN.resize(3 * nb_base_fun_row, nb_base_fun_col, false);
1769 NN.clear();
1770
1771 auto const_unit_n =
1772 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr.get()));
1773
1774 auto t_w = getFTensor0IntegrationWeightMaster();
1775
1776 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1777
1778 double val_m = t_w * area_master;
1779 auto t_base_master = row_data.getFTensor0N(gg, 0);
1780
1781 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
1782 auto t_assemble_m = get_tensor_from_mat(NN, 3 * bbr, 0);
1783 auto t_base_lambda = col_data.getFTensor0N(gg, 0);
1784 const double m = val_m * t_base_master;
1785 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
1786 const double n = m * t_base_lambda;
1787 t_assemble_m(i) -= n * const_unit_n(i);
1788 ++t_assemble_m;
1789 ++t_base_lambda; // update cols slave
1790 }
1791 ++t_base_master; // update rows master
1792 }
1793 ++t_w;
1794 }
1795
1796 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1797 ADD_VALUES);
1798 }
1799
1801}
1802
1804 int side, EntityType type, EntitiesFieldData::EntData &data) {
1806
1807 if (data.getIndices().size() == 0)
1809
1810 const int nb_dofs = data.getIndices().size();
1811 if (nb_dofs) {
1812 const int nb_gauss_pts = data.getN().size1();
1813 int nb_base_fun_col = data.getFieldData().size() / 3;
1814
1815 vecF.resize(nb_dofs, false);
1816 vecF.clear();
1817
1818 const double area_s =
1819 commonDataSimpleContact->areaSlave; // same area in master and slave
1820
1821 auto get_tensor_vec = [](VectorDouble &n, const int r) {
1822 return FTensor::Tensor1<double *, 3>(&n(r + 0), &n(r + 1), &n(r + 2));
1823 };
1824
1826
1827 auto const_unit_n = get_tensor_vec(
1828 commonDataSimpleContact->normalVectorSlavePtr.get()[0], 0);
1829
1830 auto t_w = getFTensor0IntegrationWeightSlave();
1831 auto t_aug_lambda_ptr =
1832 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1833
1834 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1835 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1836 ++t_aug_lambda_ptr;
1837 ++t_w;
1838 continue;
1839 }
1840 const double val_s = t_w * area_s;
1841
1842 FTensor::Tensor0<double *> t_base_slave(&data.getN()(gg, 0));
1843
1844 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
1845
1846 auto t_assemble_s = get_tensor_vec(vecF, 3 * bbc);
1847
1848 t_assemble_s(i) -=
1849 val_s * const_unit_n(i) * t_aug_lambda_ptr * t_base_slave;
1850
1851 ++t_base_slave;
1852 }
1853 ++t_aug_lambda_ptr;
1854 ++t_w;
1855 } // for gauss points
1856
1857 CHKERR VecSetValues(getSNESf(), data, &*vecF.begin(), ADD_VALUES);
1858 }
1860}
1861
1863 int side, EntityType type, EntitiesFieldData::EntData &data) {
1865
1866 if (data.getIndices().size() == 0)
1868
1869 const int nb_dofs = data.getIndices().size();
1870 if (nb_dofs) {
1871
1872 const int nb_gauss_pts = data.getN().size1();
1873 const int nb_base_fun_col = data.getFieldData().size() / 3;
1874
1875 vecF.resize(nb_dofs,
1876 false); // the last false in ublas
1877 // resize will destroy (not
1878 // preserved) the old
1879 // values
1880 vecF.clear();
1881
1882 const double area_m =
1883 commonDataSimpleContact->areaSlave; // same area in master and slave
1884
1885 auto get_tensor_vec = [](VectorDouble &n, const int r) {
1886 return FTensor::Tensor1<double *, 3>(&n(r + 0), &n(r + 1), &n(r + 2));
1887 };
1888
1890
1891 auto const_unit_n = get_tensor_vec(
1892 commonDataSimpleContact->normalVectorSlavePtr.get()[0], 0);
1893
1894 auto t_w = getFTensor0IntegrationWeightSlave();
1895 auto t_aug_lambda_ptr =
1896 getFTensor0FromVec(*commonDataSimpleContact->augmentedLambdasPtr);
1897
1898 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1899 if (t_aug_lambda_ptr > 0 && std::abs(t_aug_lambda_ptr) > ALM_TOL) {
1900 ++t_aug_lambda_ptr;
1901 ++t_w;
1902 continue;
1903 }
1904 const double val_m = t_w * area_m;
1905
1906 FTensor::Tensor0<double *> t_base_master(&data.getN()(gg, 0));
1907
1908 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
1909
1910 auto t_assemble_m = get_tensor_vec(vecF, 3 * bbc);
1911
1912 t_assemble_m(i) +=
1913 val_m * const_unit_n(i) * t_aug_lambda_ptr * t_base_master;
1914
1915 ++t_base_master;
1916 }
1917 ++t_aug_lambda_ptr;
1918 ++t_w;
1919 } // for gauss points
1920
1921 CHKERR VecSetValues(getSNESf(), data, &*vecF.begin(), ADD_VALUES);
1922 }
1924}
1925
1928 int row_side, int col_side, EntityType row_type, EntityType col_type,
1929 EntData &row_data, EntData &col_data) {
1931
1932 // Both sides are needed since both sides contribute their shape
1933 // function to the stiffness matrix
1934 const int nb_row = row_data.getIndices().size();
1935 const int nb_col = col_data.getIndices().size();
1936
1937 if (nb_row && nb_col) {
1938 const int nb_gauss_pts = row_data.getN().size1();
1939
1940 int nb_base_fun_row = row_data.getFieldData().size() / 3;
1941 int nb_base_fun_col = col_data.getFieldData().size();
1942
1943 const double area_slave =
1944 commonDataSimpleContact->areaSlave; // same area in master and slave
1945
1946 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
1947 return FTensor::Tensor1<double *, 3>(&m(r + 0, c), &m(r + 1, c),
1948 &m(r + 2, c));
1949 };
1950
1951 auto get_tensor_vec = [](VectorDouble &n) {
1952 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
1953 };
1954
1956
1957 NN.resize(3 * nb_base_fun_row, nb_base_fun_col, false);
1958 NN.clear();
1959
1960 auto const_unit_n =
1961 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr.get()));
1962
1963 auto t_w = getFTensor0IntegrationWeightSlave();
1964
1965 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
1966
1967 double val_m = t_w * area_slave;
1968 auto t_base_master = row_data.getFTensor0N(gg, 0);
1969
1970 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
1971 auto t_assemble_m = get_tensor_from_mat(NN, 3 * bbr, 0);
1972 auto t_base_lambda = col_data.getFTensor0N(gg, 0);
1973 const double m = val_m * t_base_master;
1974 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
1975 const double n = m * t_base_lambda;
1976 t_assemble_m(i) += n * const_unit_n(i);
1977 ++t_assemble_m;
1978 ++t_base_lambda; // update cols slave
1979 }
1980 ++t_base_master; // update rows master
1981 }
1982 ++t_w;
1983 }
1984
1985 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
1986 ADD_VALUES);
1987 }
1989}
1990
1993 int row_side, int col_side, EntityType row_type, EntityType col_type,
1994 EntData &row_data, EntData &col_data) {
1996
1997 const int nb_row = row_data.getIndices().size();
1998 const int nb_col = col_data.getIndices().size();
1999
2000 if (nb_row && nb_col) {
2001 const int nb_gauss_pts = row_data.getN().size1();
2002 const double area_slave =
2003 commonDataSimpleContact->areaSlave; // same area in master and slave
2004
2005 NN.resize(nb_row, nb_col, false);
2006 NN.clear();
2007
2008 auto t_lagrange_slave =
2009 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
2010 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
2011 auto t_w = getFTensor0IntegrationWeightSlave();
2012 const double cn_value = *cNPtr.get();
2013
2014 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2015 const double val_s = SimpleContactProblem::ConstrainFunction_dl(
2016 cn_value, t_gap_gp, t_lagrange_slave) *
2017 t_w * area_slave;
2018
2020 &*NN.data().begin());
2021
2022 auto t_base_lambda_row = row_data.getFTensor0N(gg, 0);
2023 for (int bbr = 0; bbr != nb_row; ++bbr) {
2024 auto t_base_lambda_col = col_data.getFTensor0N(gg, 0);
2025 const double s = val_s * t_base_lambda_row;
2026 for (int bbc = 0; bbc != nb_col; ++bbc) {
2027
2028 t_mat += s * t_base_lambda_col;
2029
2030 ++t_mat;
2031 ++t_base_lambda_col; // update cols
2032 }
2033 ++t_base_lambda_row; // update rows
2034 }
2035
2036 ++t_lagrange_slave;
2037 ++t_gap_gp;
2038 ++t_w;
2039 }
2040
2041 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
2042 ADD_VALUES);
2043 }
2044
2046}
2047
2050 int row_side, int col_side, EntityType row_type, EntityType col_type,
2051 EntData &row_data, EntData &col_data) {
2053
2054 const int nb_row = row_data.getIndices().size();
2055 const int nb_col = col_data.getIndices().size();
2056
2057 if (nb_row && nb_col) {
2058
2059 const int nb_gauss_pts = row_data.getN().size1();
2060 int nb_base_fun_row = row_data.getFieldData().size();
2061 int nb_base_fun_col = col_data.getFieldData().size() / 3;
2062
2063 const double area_master =
2064 commonDataSimpleContact->areaSlave; // same area in master and slave
2065
2066 const double area_slave =
2067 commonDataSimpleContact->areaSlave; // same area in master and slave
2068
2069 NN.resize(nb_base_fun_row, 3 * nb_base_fun_col, false);
2070 NN.clear();
2071
2072 auto get_tensor_from_vec = [](VectorDouble &n) {
2073 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
2074 };
2075
2077
2078 auto t_const_unit_n =
2079 get_tensor_from_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
2080
2081 auto t_const_unit_master =
2082 get_tensor_from_vec(*(commonDataSimpleContact->normalVectorMasterPtr));
2083
2084 auto t_lagrange_slave =
2085 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
2086 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
2087
2088 auto t_w = getFTensor0IntegrationWeightSlave();
2089 const double cn_value = *cNPtr.get();
2090
2091 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2092 const double val_m = SimpleContactProblem::ConstrainFunction_dg(
2093 cn_value, t_gap_gp, t_lagrange_slave) *
2094 t_w * area_slave;
2095
2096 auto t_base_lambda = row_data.getFTensor0N(gg, 0);
2097 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
2098
2099 auto t_base_master = col_data.getFTensor0N(gg, 0);
2101 &NN(bbr, 0), &NN(bbr, 1), &NN(bbr, 2)};
2102 const double m = val_m * t_base_lambda;
2103
2104 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
2105
2106 t_mat(i) += t_const_unit_n(i) * m * t_base_master;
2107
2108 ++t_base_master; // update rows
2109 ++t_mat;
2110 }
2111 ++t_base_lambda; // update cols master
2112 }
2113
2114 ++t_gap_gp;
2115 ++t_lagrange_slave;
2116 ++t_w;
2117 }
2118
2119 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
2120 ADD_VALUES);
2121 }
2122
2124}
2125
2128 int row_side, int col_side, EntityType row_type, EntityType col_type,
2129 EntData &row_data, EntData &col_data) {
2131
2132 const int nb_row = row_data.getIndices().size();
2133 const int nb_col = col_data.getIndices().size();
2134
2135 if (nb_row && nb_col) {
2136
2137 const int nb_gauss_pts = row_data.getN().size1();
2138 int nb_base_fun_row = row_data.getFieldData().size();
2139 int nb_base_fun_col = col_data.getFieldData().size() / 3;
2140
2141 const double area_slave =
2142 commonDataSimpleContact->areaSlave; // same area in master and slave
2143
2144 NN.resize(nb_base_fun_row, 3 * nb_base_fun_col, false);
2145 NN.clear();
2146
2147 auto get_tensor_from_vec = [](VectorDouble &n) {
2148 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
2149 };
2150
2152
2153 auto t_lagrange_slave =
2154 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
2155 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
2156
2157 auto t_const_unit_n =
2158 get_tensor_from_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
2159
2160 auto t_w = getFTensor0IntegrationWeightSlave();
2161 const double cn_value = *cNPtr.get();
2162 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2163 const double val_m = SimpleContactProblem::ConstrainFunction_dg(
2164 cn_value, t_gap_gp, t_lagrange_slave) *
2165 t_w * area_slave;
2166
2167 auto t_base_lambda = row_data.getFTensor0N(gg, 0);
2168 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
2169
2170 auto t_base_slave = col_data.getFTensor0N(gg, 0);
2172 &NN(bbr, 0), &NN(bbr, 1), &NN(bbr, 2)};
2173 const double m = val_m * t_base_lambda;
2174
2175 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
2176
2177 t_mat(i) -= t_const_unit_n(i) * m * t_base_slave;
2178
2179 ++t_base_slave; // update rows
2180 ++t_mat;
2181 }
2182 ++t_base_lambda; // update cols master
2183 }
2184
2185 ++t_gap_gp;
2186 ++t_lagrange_slave;
2187 ++t_w;
2188 }
2189
2190 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*NN.data().begin(),
2191 ADD_VALUES);
2192 }
2193
2195}
2196
2198 EntityType type,
2199 EntData &data) {
2201
2202 if (type != MBVERTEX)
2204
2205 const EntityHandle prism_ent = getFEEntityHandle();
2206 EntityHandle tri_ent;
2207 if (stateTagSide == MASTER_SIDE) {
2208 CHKERR mField.get_moab().side_element(prism_ent, 2, 3, tri_ent);
2209 }
2210 if (stateTagSide == SLAVE_SIDE) {
2211 CHKERR mField.get_moab().side_element(prism_ent, 2, 4, tri_ent);
2212 }
2213
2214 int nb_dofs = data.getFieldData().size();
2215 if (nb_dofs == 0)
2217 int nb_gauss_pts = data.getN().size1();
2218
2219 double def_val = 0.;
2220
2221 Tag th_gap;
2222 CHKERR moabOut.tag_get_handle("GAP", 1, MB_TYPE_DOUBLE, th_gap,
2223 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
2224
2225 Tag th_lag_mult;
2226 CHKERR moabOut.tag_get_handle("LAGMULT", 1, MB_TYPE_DOUBLE, th_lag_mult,
2227 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
2228
2229 Tag th_lag_gap_prod;
2230 CHKERR moabOut.tag_get_handle("LAG_GAP_PROD", 1, MB_TYPE_DOUBLE,
2231 th_lag_gap_prod, MB_TAG_CREAT | MB_TAG_SPARSE,
2232 &def_val);
2233
2234 int def_val_int = 0;
2235
2236 Tag th_state;
2237 CHKERR moabOut.tag_get_handle("STATE", 1, MB_TYPE_INTEGER, th_state,
2238 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
2239
2240 Tag th_state_side;
2241 if (stateTagSide > 0) {
2242 CHKERR mField.get_moab().tag_get_handle(
2243 "STATE", 1, MB_TYPE_INTEGER, th_state_side,
2244 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
2245 }
2246
2247 auto get_tag_pos = [&](const std::string name) {
2248 Tag th;
2249 constexpr std::array<double, 3> def_vals = {0, 0, 0};
2250 CHKERR moabOut.tag_get_handle(name.c_str(), 3, MB_TYPE_DOUBLE, th,
2251 MB_TAG_CREAT | MB_TAG_SPARSE,
2252 def_vals.data());
2253 return th;
2254 };
2255 auto th_pos_master = get_tag_pos("MASTER_SPATIAL_POSITION");
2256 auto th_pos_slave = get_tag_pos("SLAVE_SPATIAL_POSITION");
2257 auto th_master_coords = get_tag_pos("MASTER_GAUSS_PTS_COORDS");
2258
2259 // FIXME: why getFEEntityHandle ?
2260 EntityHandle new_vertex = getFEEntityHandle();
2261
2262 auto t_gap_ptr = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
2263
2264 auto t_lagrange_slave =
2265 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
2266 auto t_lag_gap_prod_slave =
2267 getFTensor0FromVec(*commonDataSimpleContact->lagGapProdPtr);
2268 auto t_position_master = getFTensor1FromMat<3>(
2269 *commonDataSimpleContact->positionAtGaussPtsMasterPtr);
2270 auto t_position_slave = getFTensor1FromMat<3>(
2271 *commonDataSimpleContact->positionAtGaussPtsSlavePtr);
2272
2273 auto get_ftensor_coords_at_gauss_pts_slave = [&](auto &coords_at_gauss_pts) {
2274 auto ptr = &*coords_at_gauss_pts.data().begin();
2275 return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(&ptr[0], &ptr[1],
2276 &ptr[2]);
2277 };
2278 auto t_coords_at_gauss_pts_slave =
2279 get_ftensor_coords_at_gauss_pts_slave(getCoordsAtGaussPtsSlave());
2280 auto t_coords_at_gauss_pts_master =
2281 get_ftensor_coords_at_gauss_pts_slave(getCoordsAtGaussPtsMaster());
2282
2283 auto t_state_ptr =
2284 getFTensor0FromVec(*commonDataSimpleContact->gaussPtsStatePtr);
2285
2286 auto set_float_precision = [](const double x) {
2287 if (std::abs(x) < std::numeric_limits<float>::epsilon())
2288 return 0.;
2289 else
2290 return x;
2291 };
2292
2293 std::array<double, 3> pos_vec;
2294
2295 auto get_vec_ptr = [&](auto t) {
2296 for (int dd = 0; dd != 3; ++dd)
2297 pos_vec[dd] = set_float_precision(t(dd));
2298 return pos_vec.data();
2299 };
2300
2301 int count_active_pts = 0;
2302
2303 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2304
2305 int state = 0;
2306 if (t_state_ptr > 0.5) {
2307 state = 1;
2308 ++count_active_pts;
2309 }
2310
2311 bool output_at_gauss_pts = true;
2312 if (!postProcSurface.empty()) {
2313 Range tri_ents;
2314 CHKERR mField.get_moab().get_adjacencies(
2315 &prism_ent, 1, 2, false, tri_ents, moab::Interface::UNION);
2316 tri_ents = tri_ents.subset_by_type(MBTRI);
2317 if (intersect(postProcSurface, tri_ents).empty())
2318 output_at_gauss_pts = false;
2319 }
2320
2321 if (output_at_gauss_pts) {
2322 CHKERR moabOut.create_vertex(get_vec_ptr(t_coords_at_gauss_pts_slave),
2323 new_vertex);
2324
2325 double gap_vtk = set_float_precision(t_gap_ptr);
2326 CHKERR moabOut.tag_set_data(th_gap, &new_vertex, 1, &gap_vtk);
2327
2328 double lag_gap_prod_vtk = set_float_precision(t_lag_gap_prod_slave);
2329 CHKERR moabOut.tag_set_data(th_lag_gap_prod, &new_vertex, 1,
2330 &lag_gap_prod_vtk);
2331
2332 double lagrange_slave_vtk = set_float_precision(t_lagrange_slave);
2333 CHKERR moabOut.tag_set_data(th_lag_mult, &new_vertex, 1,
2334 &lagrange_slave_vtk);
2335
2336 CHKERR moabOut.tag_set_data(th_state, &new_vertex, 1, &state);
2337
2338 CHKERR moabOut.tag_set_data(th_pos_master, &new_vertex, 1,
2339 get_vec_ptr(t_position_master));
2340 CHKERR moabOut.tag_set_data(th_pos_slave, &new_vertex, 1,
2341 get_vec_ptr(t_position_slave));
2342
2343 CHKERR moabOut.tag_set_data(th_master_coords, &new_vertex, 1,
2344 get_vec_ptr(t_coords_at_gauss_pts_master));
2345 }
2346
2347 ++t_gap_ptr;
2348 ++t_lagrange_slave;
2349 ++t_lag_gap_prod_slave;
2350 ++t_position_master;
2351 ++t_position_slave;
2352 ++t_coords_at_gauss_pts_slave;
2353 ++t_coords_at_gauss_pts_master;
2354 ++t_state_ptr;
2355 }
2356
2357 if (stateTagSide > 0) {
2358 int state_side = 0;
2359 if (count_active_pts >= nb_gauss_pts / 2) {
2360 state_side = 1;
2361 }
2362 CHKERR mField.get_moab().tag_set_data(th_state_side, &tri_ent, 1,
2363 &state_side);
2364 }
2365
2367}
2368
2370 EntityType type,
2371 EntData &data) {
2373
2374 if (type != MBVERTEX)
2376
2377 int nb_dofs = data.getFieldData().size();
2378 if (nb_dofs == 0)
2380 int nb_gauss_pts = data.getN().size1();
2381
2382 auto t_gap_ptr = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
2383
2384 auto t_lagrange_slave =
2385 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
2386 double d_gap;
2387 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2388 const double d_lambda =
2389 std::abs(t_lagrange_slave) < TOL ? 0.0 : t_lagrange_slave;
2390 d_gap = std::abs(t_gap_ptr) < TOL ? 0.0 : t_gap_ptr;
2391 mySplit << d_lambda << " " << d_gap << " " << std::endl;
2392
2393 ++t_gap_ptr;
2394 ++t_lagrange_slave;
2395 }
2397}
2398
2400 boost::shared_ptr<SimpleContactElement> fe_rhs_simple_contact,
2401 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
2402 string field_name, string lagrange_field_name, bool is_alm,
2403 bool is_eigen_pos_field, string eigen_pos_field_name,
2404 bool use_reference_coordinates) {
2406
2407 fe_rhs_simple_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
2408 "MESH_NODE_POSITIONS", common_data_simple_contact));
2409
2410 fe_rhs_simple_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
2411 "MESH_NODE_POSITIONS", common_data_simple_contact));
2412
2413 fe_rhs_simple_contact->getOpPtrVector().push_back(
2415 common_data_simple_contact));
2416
2417 fe_rhs_simple_contact->getOpPtrVector().push_back(
2418 new OpGetPositionAtGaussPtsSlave(field_name, common_data_simple_contact));
2419
2420 if (is_eigen_pos_field) {
2421 fe_rhs_simple_contact->getOpPtrVector().push_back(
2423 eigen_pos_field_name, common_data_simple_contact));
2424
2425 fe_rhs_simple_contact->getOpPtrVector().push_back(
2427 eigen_pos_field_name, common_data_simple_contact));
2428 if (use_reference_coordinates) {
2429 fe_rhs_simple_contact->getOpPtrVector().push_back(
2430 new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
2431 common_data_simple_contact));
2432
2433 fe_rhs_simple_contact->getOpPtrVector().push_back(
2434 new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
2435 common_data_simple_contact));
2436 }
2437 }
2438
2439 fe_rhs_simple_contact->getOpPtrVector().push_back(
2440 new OpGetGapSlave(field_name, common_data_simple_contact));
2441
2442 fe_rhs_simple_contact->getOpPtrVector().push_back(
2443 new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
2444 common_data_simple_contact));
2445
2446 fe_rhs_simple_contact->getOpPtrVector().push_back(new OpGetGaussPtsState(
2447 lagrange_field_name, common_data_simple_contact, cnValue, is_alm));
2448
2449 if (!is_alm) {
2450 fe_rhs_simple_contact->getOpPtrVector().push_back(
2452 common_data_simple_contact));
2453
2454 fe_rhs_simple_contact->getOpPtrVector().push_back(new OpCalIntCompFunSlave(
2455 lagrange_field_name, common_data_simple_contact, cnValuePtr));
2456 } else {
2457 fe_rhs_simple_contact->getOpPtrVector().push_back(
2458 new OpGetAugmentedLambdaSlave(field_name, common_data_simple_contact,
2459 cnValue));
2460
2461 fe_rhs_simple_contact->getOpPtrVector().push_back(
2463 common_data_simple_contact));
2464
2465 fe_rhs_simple_contact->getOpPtrVector().push_back(
2466 new OpGapConstraintAugmentedRhs(lagrange_field_name,
2467 common_data_simple_contact, cnValue));
2468 }
2469
2471}
2472
2474 boost::shared_ptr<SimpleContactElement> fe_rhs_simple_contact,
2475 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
2476 string field_name, string lagrange_field_name, bool is_alm,
2477 bool is_eigen_pos_field, string eigen_pos_field_name,
2478 bool use_reference_coordinates) {
2480
2481 fe_rhs_simple_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
2482 "MESH_NODE_POSITIONS", common_data_simple_contact));
2483
2484 fe_rhs_simple_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
2485 "MESH_NODE_POSITIONS", common_data_simple_contact));
2486
2487 fe_rhs_simple_contact->getOpPtrVector().push_back(
2488 new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
2489 common_data_simple_contact));
2490
2491 if (!is_alm) {
2492 fe_rhs_simple_contact->getOpPtrVector().push_back(
2494 common_data_simple_contact));
2495 } else {
2496
2497 fe_rhs_simple_contact->getOpPtrVector().push_back(
2499 common_data_simple_contact));
2500
2501 fe_rhs_simple_contact->getOpPtrVector().push_back(
2503 common_data_simple_contact));
2504
2505 if (is_eigen_pos_field) {
2506 fe_rhs_simple_contact->getOpPtrVector().push_back(
2508 eigen_pos_field_name, common_data_simple_contact));
2509
2510 fe_rhs_simple_contact->getOpPtrVector().push_back(
2512 eigen_pos_field_name, common_data_simple_contact));
2513
2514 if (use_reference_coordinates) {
2515 fe_rhs_simple_contact->getOpPtrVector().push_back(
2517 "MESH_NODE_POSITIONS", common_data_simple_contact));
2518
2519 fe_rhs_simple_contact->getOpPtrVector().push_back(
2520 new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
2521 common_data_simple_contact));
2522 }
2523 }
2524
2525 fe_rhs_simple_contact->getOpPtrVector().push_back(
2526 new OpGetGapSlave(field_name, common_data_simple_contact));
2527
2528 fe_rhs_simple_contact->getOpPtrVector().push_back(
2529 new OpGetAugmentedLambdaSlave(field_name, common_data_simple_contact,
2530 cnValue));
2531
2532 fe_rhs_simple_contact->getOpPtrVector().push_back(
2534 common_data_simple_contact));
2535 }
2536
2538}
2539
2541 boost::shared_ptr<SimpleContactElement> fe_lhs_simple_contact,
2542 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
2543 string field_name, string lagrange_field_name, bool is_alm,
2544 bool is_eigen_pos_field, string eigen_pos_field_name,
2545 bool use_reference_coordinates) {
2547
2548 fe_lhs_simple_contact->getOpPtrVector().push_back(
2550 common_data_simple_contact));
2551
2552 fe_lhs_simple_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
2553 "MESH_NODE_POSITIONS", common_data_simple_contact));
2554
2555 fe_lhs_simple_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
2556 "MESH_NODE_POSITIONS", common_data_simple_contact));
2557
2558 fe_lhs_simple_contact->getOpPtrVector().push_back(
2559 new OpGetPositionAtGaussPtsSlave(field_name, common_data_simple_contact));
2560
2561 if (is_eigen_pos_field) {
2562 fe_lhs_simple_contact->getOpPtrVector().push_back(
2564 eigen_pos_field_name, common_data_simple_contact));
2565
2566 fe_lhs_simple_contact->getOpPtrVector().push_back(
2568 eigen_pos_field_name, common_data_simple_contact));
2569
2570 if (use_reference_coordinates) {
2571 fe_lhs_simple_contact->getOpPtrVector().push_back(
2572 new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
2573 common_data_simple_contact));
2574
2575 fe_lhs_simple_contact->getOpPtrVector().push_back(
2576 new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
2577 common_data_simple_contact));
2578 }
2579 }
2580
2581 fe_lhs_simple_contact->getOpPtrVector().push_back(
2582 new OpGetGapSlave(field_name, common_data_simple_contact));
2583
2584 fe_lhs_simple_contact->getOpPtrVector().push_back(
2585 new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
2586 common_data_simple_contact));
2587 if (!is_alm) {
2588 fe_lhs_simple_contact->getOpPtrVector().push_back(
2590 field_name, lagrange_field_name, common_data_simple_contact));
2591
2592 fe_lhs_simple_contact->getOpPtrVector().push_back(
2594 lagrange_field_name, common_data_simple_contact, cnValuePtr));
2595
2596 fe_lhs_simple_contact->getOpPtrVector().push_back(
2598 lagrange_field_name, field_name, common_data_simple_contact,
2599 cnValuePtr));
2600
2601 fe_lhs_simple_contact->getOpPtrVector().push_back(
2603 lagrange_field_name, field_name, common_data_simple_contact,
2604 cnValuePtr));
2605 } else {
2606 fe_lhs_simple_contact->getOpPtrVector().push_back(
2607 new OpGetAugmentedLambdaSlave(field_name, common_data_simple_contact,
2608 cnValue));
2609
2610 fe_lhs_simple_contact->getOpPtrVector().push_back(
2612 field_name, lagrange_field_name, common_data_simple_contact));
2613
2614 fe_lhs_simple_contact->getOpPtrVector().push_back(
2616 field_name, field_name, cnValue, common_data_simple_contact));
2617
2618 fe_lhs_simple_contact->getOpPtrVector().push_back(
2620 field_name, field_name, cnValue, common_data_simple_contact));
2621
2622 fe_lhs_simple_contact->getOpPtrVector().push_back(
2624 lagrange_field_name, common_data_simple_contact, cnValue));
2625
2626 fe_lhs_simple_contact->getOpPtrVector().push_back(
2628 field_name, lagrange_field_name, common_data_simple_contact,
2629 cnValue));
2630
2631 fe_lhs_simple_contact->getOpPtrVector().push_back(
2633 field_name, lagrange_field_name, common_data_simple_contact,
2634 cnValue));
2635 }
2636
2638}
2639
2641 boost::shared_ptr<SimpleContactElement> fe_lhs_simple_contact,
2642 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
2643 string field_name, string lagrange_field_name, bool is_alm,
2644 bool is_eigen_pos_field, string eigen_pos_field_name,
2645 bool use_reference_coordinates) {
2647
2648 fe_lhs_simple_contact->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
2649 "MESH_NODE_POSITIONS", common_data_simple_contact));
2650
2651 fe_lhs_simple_contact->getOpPtrVector().push_back(new OpGetNormalMasterALE(
2652 "MESH_NODE_POSITIONS", common_data_simple_contact));
2653
2654 fe_lhs_simple_contact->getOpPtrVector().push_back(
2655 new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
2656 common_data_simple_contact));
2657 if (!is_alm) {
2658 fe_lhs_simple_contact->getOpPtrVector().push_back(
2660 field_name, lagrange_field_name, common_data_simple_contact));
2661 } else {
2662
2663 fe_lhs_simple_contact->getOpPtrVector().push_back(
2665 common_data_simple_contact));
2666
2667 fe_lhs_simple_contact->getOpPtrVector().push_back(
2669 common_data_simple_contact));
2670
2671 if (is_eigen_pos_field) {
2672 fe_lhs_simple_contact->getOpPtrVector().push_back(
2674 eigen_pos_field_name, common_data_simple_contact));
2675
2676 fe_lhs_simple_contact->getOpPtrVector().push_back(
2678 eigen_pos_field_name, common_data_simple_contact));
2679 if (use_reference_coordinates) {
2680 fe_lhs_simple_contact->getOpPtrVector().push_back(
2682 "MESH_NODE_POSITIONS", common_data_simple_contact));
2683
2684 fe_lhs_simple_contact->getOpPtrVector().push_back(
2685 new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
2686 common_data_simple_contact));
2687 }
2688 }
2689
2690 fe_lhs_simple_contact->getOpPtrVector().push_back(
2691 new OpGetGapSlave(field_name, common_data_simple_contact));
2692
2693 fe_lhs_simple_contact->getOpPtrVector().push_back(
2694 new OpGetAugmentedLambdaSlave(field_name, common_data_simple_contact,
2695 cnValue));
2696
2697 fe_lhs_simple_contact->getOpPtrVector().push_back(
2699 field_name, lagrange_field_name, common_data_simple_contact));
2700
2701 fe_lhs_simple_contact->getOpPtrVector().push_back(
2703 field_name, field_name, cnValue, common_data_simple_contact));
2704
2705 fe_lhs_simple_contact->getOpPtrVector().push_back(
2707 field_name, field_name, cnValue, common_data_simple_contact));
2708 }
2709
2711}
2712
2714 boost::shared_ptr<ConvectMasterContactElement> fe_lhs_simple_contact,
2715 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
2716 string field_name, string lagrange_field_name, bool is_alm,
2717 bool is_eigen_pos_field, string eigen_pos_field_name,
2718 bool use_reference_coordinates) {
2721 boost::dynamic_pointer_cast<SimpleContactElement>(fe_lhs_simple_contact),
2722 common_data_simple_contact, field_name, lagrange_field_name, is_alm,
2723 is_eigen_pos_field, eigen_pos_field_name);
2724
2725 fe_lhs_simple_contact->getOpPtrVector().push_back(
2726 new OpCalculateGradPositionXi(field_name, common_data_simple_contact));
2727
2728 fe_lhs_simple_contact->getOpPtrVector().push_back(
2730 lagrange_field_name, field_name, common_data_simple_contact,
2731 cnValuePtr, ContactOp::FACESLAVESLAVE,
2732 fe_lhs_simple_contact->getConvectPtr()->getDiffKsiSpatialSlave()));
2733
2734 fe_lhs_simple_contact->getOpPtrVector().push_back(
2736 lagrange_field_name, field_name, common_data_simple_contact,
2737 cnValuePtr, ContactOp::FACESLAVEMASTER,
2738 fe_lhs_simple_contact->getConvectPtr()->getDiffKsiSpatialMaster()));
2739
2741}
2742
2744 boost::shared_ptr<ConvectSlaveContactElement> fe_lhs_simple_contact,
2745 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
2746 string field_name, string lagrange_field_name, bool is_alm,
2747 bool is_eigen_pos_field, string eigen_pos_field_name,
2748 bool use_reference_coordinates) {
2750
2752 boost::dynamic_pointer_cast<SimpleContactElement>(fe_lhs_simple_contact),
2753 common_data_simple_contact, field_name, lagrange_field_name);
2754
2755 fe_lhs_simple_contact->getOpPtrVector().push_back(new OpCalculateGradLambdaXi(
2756 lagrange_field_name, common_data_simple_contact));
2757
2758 fe_lhs_simple_contact->getOpPtrVector().push_back(
2760 field_name, field_name, common_data_simple_contact,
2761 ContactOp::FACEMASTERSLAVE,
2762 fe_lhs_simple_contact->getConvectPtr()->getDiffKsiSpatialSlave()));
2763
2764 fe_lhs_simple_contact->getOpPtrVector().push_back(
2766 field_name, field_name, common_data_simple_contact,
2767 ContactOp::FACEMASTERMASTER,
2768 fe_lhs_simple_contact->getConvectPtr()->getDiffKsiSpatialMaster()));
2770}
2771
2773 boost::shared_ptr<SimpleContactElement> fe_post_proc_simple_contact,
2774 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
2775 MoFEM::Interface &m_field, string field_name, string lagrange_field_name,
2776 moab::Interface &moab_out, bool alm_flag, bool is_eigen_pos_field,
2777 string eigen_pos_field_name, bool use_reference_coordinates,
2778 StateTagSide state_tag_side) {
2780
2781 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2782 new OpGetNormalMasterALE("MESH_NODE_POSITIONS",
2783 common_data_simple_contact));
2784
2785 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2786 new OpGetNormalSlaveALE("MESH_NODE_POSITIONS",
2787 common_data_simple_contact));
2788
2789 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2791 common_data_simple_contact));
2792
2793 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2794 new OpGetPositionAtGaussPtsSlave(field_name, common_data_simple_contact));
2795
2796 if (is_eigen_pos_field) {
2797 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2799 eigen_pos_field_name, common_data_simple_contact));
2800
2801 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2803 eigen_pos_field_name, common_data_simple_contact));
2804 if (use_reference_coordinates) {
2805 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2806 new OpGetMatPosForDisplAtGaussPtsMaster("MESH_NODE_POSITIONS",
2807 common_data_simple_contact));
2808
2809 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2810 new OpGetMatPosForDisplAtGaussPtsSlave("MESH_NODE_POSITIONS",
2811 common_data_simple_contact));
2812 }
2813 }
2814
2815 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2816 new OpGetGapSlave(field_name, common_data_simple_contact));
2817
2818 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2819 new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
2820 common_data_simple_contact));
2821
2822 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2823 new OpLagGapProdGaussPtsSlave(lagrange_field_name,
2824 common_data_simple_contact));
2825
2826 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2827 new OpGetGaussPtsState(lagrange_field_name, common_data_simple_contact,
2828 cnValue, alm_flag));
2829
2830 fe_post_proc_simple_contact->getOpPtrVector().push_back(
2831 new OpMakeVtkSlave(m_field, field_name, common_data_simple_contact,
2832 moab_out, state_tag_side));
2833
2834 fe_post_proc_simple_contact->getOpPtrVector().push_back(new OpGetContactArea(
2835 lagrange_field_name, common_data_simple_contact, cnValue, alm_flag));
2836
2838}
2839
2842 EntData &data) {
2844 const int nb_dofs = data.getFieldData().size();
2845 const int nb_integration_pts = getGaussPtsSlave().size2();
2846 auto &xi_grad_mat = *(commonDataSimpleContact->gradKsiLambdaAtGaussPtsPtr);
2847 xi_grad_mat.resize(2, nb_integration_pts, false);
2848 if (type == MBVERTEX)
2849 xi_grad_mat.clear();
2850
2852
2853 if (nb_dofs) {
2854
2855 auto t_diff_lambda_xi = getFTensor1FromMat<2>(xi_grad_mat);
2856
2857 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
2858 auto t_data = data.getFTensor0FieldData();
2859 auto t_diff_base = data.getFTensor1DiffN<2>(gg, 0);
2860 for (size_t bb = 0; bb != nb_dofs; ++bb) {
2861 t_diff_lambda_xi(I) += t_diff_base(I) * t_data;
2862 ++t_data;
2863 ++t_diff_base;
2864 }
2865 ++t_diff_lambda_xi;
2866 }
2867 }
2868
2870}
2871
2874 int row_side, int col_side, EntityType row_type, EntityType col_type,
2875 EntData &row_data, EntData &col_data) {
2877
2878 const int nb_row_dofs = row_data.getIndices().size();
2879 const int nb_col_dofs = col_data.getIndices().size();
2880 if (nb_row_dofs && (nb_col_dofs && col_type == MBVERTEX)) {
2881
2882 const int nb_gauss_pts = getGaussPtsSlave().size2();
2883 int nb_base_fun_row = nb_row_dofs / 3;
2884 int nb_base_fun_col = nb_col_dofs / 3;
2885
2886 matLhs.resize(nb_row_dofs, nb_col_dofs, false);
2887 matLhs.clear();
2888
2889 const double area_s =
2890 commonDataSimpleContact->areaSlave; // same area in master and slave
2891
2895
2896 auto get_tensor_vec = [](VectorDouble &n) {
2897 return FTensor::Tensor1<double *, 3>(&n[0], &n[1], &n[2]);
2898 };
2899
2900 auto t_const_unit_n =
2901 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
2902 auto t_diff_lambda_xi = getFTensor1FromMat<2>(
2903 *(commonDataSimpleContact->gradKsiLambdaAtGaussPtsPtr));
2904 auto t_w = getFTensor0IntegrationWeightSlave();
2905
2906 auto get_diff_ksi = [](auto &m, auto gg) {
2908 &m(0, gg), &m(1, gg), &m(2, gg), &m(3, gg), &m(4, gg), &m(5, gg));
2909 };
2910
2911 auto t_base_row = row_data.getFTensor0N();
2912
2913 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
2914
2915 double val_s = t_w * area_s;
2916 auto t_base_row = row_data.getFTensor0N(gg, 0);
2917
2918 for (int rr = 0; rr != nb_base_fun_row; ++rr) {
2919
2921 &matLhs(3 * rr + 0, 0), &matLhs(3 * rr + 0, 1),
2922 &matLhs(3 * rr + 0, 2),
2923
2924 &matLhs(3 * rr + 1, 0), &matLhs(3 * rr + 1, 1),
2925 &matLhs(3 * rr + 1, 2),
2926
2927 &matLhs(3 * rr + 2, 0), &matLhs(3 * rr + 2, 1),
2928 &matLhs(3 * rr + 2, 2)};
2929
2930 auto t_diff_convect = get_diff_ksi(*diffConvect, 3 * gg);
2931
2932 for (int cc = 0; cc != nb_base_fun_col; ++cc) {
2933 t_mat(i, j) -= val_s * t_base_row * t_const_unit_n(i) *
2934 (t_diff_lambda_xi(I) * t_diff_convect(I, j));
2935
2936 ++t_diff_convect;
2937 ++t_mat;
2938 }
2939
2940 ++t_base_row;
2941 }
2942
2943 ++t_diff_lambda_xi;
2944 ++t_w;
2945 } // for gauss points
2946
2947 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*matLhs.data().begin(),
2948 ADD_VALUES);
2949 }
2950
2952}
2953
2955 int side, EntityType type, EntData &data) {
2957 const int nb_dofs = data.getFieldData().size();
2958 const int nb_integration_pts = getGaussPtsSlave().size2();
2959 auto &xi_grad_mat = *(commonDataSimpleContact->gradKsiPositionAtGaussPtsPtr);
2960 xi_grad_mat.resize(6, nb_integration_pts, false);
2961 if (type == MBVERTEX)
2962 xi_grad_mat.clear();
2963
2966
2967 if (nb_dofs) {
2968
2969 auto t_grad_pos_xi = getFTensor2FromMat<3, 2>(xi_grad_mat);
2970
2971 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
2972 auto t_data = data.getFTensor1FieldData<3>();
2973 auto t_diff_base = data.getFTensor1DiffN<2>(gg, 0);
2974 for (size_t bb = 0; bb != nb_dofs / 3; ++bb) {
2975 t_grad_pos_xi(i, I) += t_diff_base(I) * t_data(i);
2976 ++t_data;
2977 ++t_diff_base;
2978 }
2979 ++t_grad_pos_xi;
2980 }
2981 }
2982
2984}
2985
2988 int row_side, int col_side, EntityType row_type, EntityType col_type,
2989 EntData &row_data, EntData &col_data) {
2991
2992 const int nb_row = row_data.getIndices().size();
2993 const int nb_col = col_data.getIndices().size();
2994
2995 if (nb_row && nb_col && col_type == MBVERTEX) {
2996
2997 const int nb_gauss_pts = row_data.getN().size1();
2998 int nb_base_fun_row = nb_row;
2999 int nb_base_fun_col = nb_col / 3;
3000
3001 const double area_slave =
3002 commonDataSimpleContact->areaSlave; // same area in master and slave
3003
3004 matLhs.resize(nb_row, nb_col, false);
3005 matLhs.clear();
3006
3007 auto get_diff_ksi = [](auto &m, auto gg) {
3009 &m(0, gg), &m(1, gg), &m(2, gg), &m(3, gg), &m(4, gg), &m(5, gg));
3010 };
3011
3012 auto get_tensor_from_vec = [](VectorDouble &n) {
3013 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3014 };
3015
3019
3020 auto t_const_unit_n =
3021 get_tensor_from_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
3022
3023 auto t_lagrange_slave =
3024 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
3025 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
3026 auto &xi_grad_mat =
3027 *(commonDataSimpleContact->gradKsiPositionAtGaussPtsPtr);
3028 auto t_grad = getFTensor2FromMat<3, 2>(xi_grad_mat);
3029
3030 auto t_w = getFTensor0IntegrationWeightSlave();
3031 const double cn_value = *cNPtr.get();
3032 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
3033 const double val_m = SimpleContactProblem::ConstrainFunction(
3034 cn_value, t_gap_gp, t_lagrange_slave) *
3035 t_w * area_slave;
3036 const double val_diff_m_l = SimpleContactProblem::ConstrainFunction_dl(
3037 cn_value, t_gap_gp, t_lagrange_slave) *
3038 t_w * area_slave;
3039 const double val_diff_m_g = SimpleContactProblem::ConstrainFunction_dg(
3040 cn_value, t_gap_gp, t_lagrange_slave) *
3041 t_w * area_slave;
3042
3043 auto t_base_lambda = row_data.getFTensor0N(gg, 0);
3044 auto t_base_diff_lambda = row_data.getFTensor1DiffN<2>(gg, 0);
3045
3047 &matLhs(0, 0), &matLhs(0, 1), &matLhs(0, 2)};
3048
3049 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
3050
3051 auto t_base_diff_disp = col_data.getFTensor1DiffN<2>(gg, 0);
3052 auto t_diff_convect = get_diff_ksi(*diffConvect, 3 * gg);
3053
3054 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
3055
3056 t_mat(i) += t_base_lambda * val_diff_m_g * t_const_unit_n(j) *
3057 t_grad(j, I) * t_diff_convect(I, i);
3058
3059 ++t_base_diff_disp;
3060 ++t_diff_convect;
3061 ++t_mat;
3062 }
3063
3064 ++t_base_lambda;
3065 ++t_base_diff_lambda;
3066 }
3067
3068 ++t_gap_gp;
3069 ++t_lagrange_slave;
3070 ++t_grad;
3071 ++t_w;
3072 }
3073
3074 CHKERR MatSetValues(getSNESB(), row_data, col_data, &*matLhs.data().begin(),
3075 ADD_VALUES);
3076 }
3078}
3079
3082 EntData &row_data) {
3083
3088
3089 // get number of integration points
3090 const int nb_integration_pts = getGaussPts().size2();
3091
3092 auto t_h = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->hMat);
3093 auto t_H = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->HMat);
3094
3095 commonDataSimpleContact->detHVec->resize(nb_integration_pts, false);
3096 commonDataSimpleContact->invHMat->resize(9, nb_integration_pts, false);
3097 commonDataSimpleContact->FMat->resize(9, nb_integration_pts, false);
3098
3099 commonDataSimpleContact->detHVec->clear();
3100 commonDataSimpleContact->invHMat->clear();
3101 commonDataSimpleContact->FMat->clear();
3102
3103 auto t_detH = getFTensor0FromVec(*commonDataSimpleContact->detHVec);
3104 auto t_invH = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->invHMat);
3105 auto t_F = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->FMat);
3106
3107 for (int gg = 0; gg != nb_integration_pts; ++gg) {
3108 CHKERR determinantTensor3by3(t_H, t_detH);
3109 CHKERR invertTensor3by3(t_H, t_detH, t_invH);
3110 t_F(i, j) = t_h(i, k) * t_invH(k, j);
3111 ++t_h;
3112 ++t_H;
3113 ++t_detH;
3114 ++t_invH;
3115 ++t_F;
3116 }
3117
3119}
3120
3122 int side, EntityType type, EntData &data) {
3124
3125 if (type != MBTRI)
3127
3128 int side_number;
3129 if (faceType == ContactOp::FACEMASTER)
3130 side_number = 3;
3131 else
3132 side_number = 4;
3133
3134 const EntityHandle tri = getSideEntity(side_number, type);
3135
3136 CHKERR loopSideVolumes(sideFeName, *sideFe, 3, tri);
3137
3139}
3140
3143 EntData &row_data) {
3144
3146
3147 // get number of dofs on row
3148 nbRows = row_data.getIndices().size();
3149 // if no dofs on row, exit that work, nothing to do here
3150 if (!nbRows)
3152
3153 vecF.resize(nbRows, false);
3154 vecF.clear();
3155
3156 // get number of integration points
3157 nbIntegrationPts = getGaussPtsMaster().size2();
3158
3159 // integrate local matrix for entity block
3160 CHKERR iNtegrate(row_data);
3161
3162 // assemble local matrix
3163 CHKERR aSsemble(row_data);
3164
3166}
3167
3171
3172 if (data.getIndices().size() == 0)
3174 // in case the tet is not in database
3175 if (commonDataSimpleContact->FMat->size1() != 9)
3177
3178 const int nb_base_fun_col = nbRows / 3;
3179
3180 auto get_tensor_vec = [](VectorDouble &n, const int r) {
3181 return FTensor::Tensor1<double *, 3>(&n(r + 0), &n(r + 1), &n(r + 2));
3182 };
3183
3186 auto lagrange_slave =
3187 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
3188
3189 auto t_F = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->FMat);
3190 auto normal_at_gp =
3191 get_tensor_vec(*(commonDataSimpleContact->normalVectorMasterPtr), 0);
3192
3193 auto t_w = getFTensor0IntegrationWeightMaster();
3194 double &area_master = commonDataSimpleContact->areaMaster;
3195 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
3196 const double val_m = area_master * t_w;
3197 FTensor::Tensor0<double *> t_base_master(&data.getN()(gg, 0));
3198
3199 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
3200
3201 const double s = val_m * t_base_master * lagrange_slave;
3202
3203 auto t_assemble_s = get_tensor_vec(vecF, 3 * bbc);
3204
3205 t_assemble_s(i) -= s * t_F(j, i) * normal_at_gp(j);
3206
3207 ++t_base_master;
3208 }
3209 ++t_F;
3210 ++lagrange_slave;
3211 ++t_w;
3212 } // for gauss points
3213
3215}
3216
3220
3221 // get pointer to first global index on row
3222 const int *row_indices = &*row_data.getIndices().data().begin();
3223 auto &data = *commonDataSimpleContact;
3224 if (data.forcesOnlyOnEntitiesRow.empty())
3226
3227 rowIndices.resize(nbRows, false);
3228 noalias(rowIndices) = row_data.getIndices();
3229 row_indices = &rowIndices[0];
3230 VectorDofs &dofs = row_data.getFieldDofs();
3231 VectorDofs::iterator dit = dofs.begin();
3232 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
3233 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
3234 data.forcesOnlyOnEntitiesRow.end()) {
3235 rowIndices[ii] = -1;
3236 }
3237 }
3238
3239 CHKERR VecSetValues(getSNESf(), nbRows, row_indices, &*vecF.data().begin(),
3240 ADD_VALUES);
3242}
3243
3245SimpleContactProblem::OpCalMatForcesALESlave::doWork(int side, EntityType type,
3246 EntData &row_data) {
3247
3249
3250 // get number of dofs on row
3251 nbRows = row_data.getIndices().size();
3252 // if no dofs on row, exit that work, nothing to do here
3253 if (!nbRows)
3255
3256 vecF.resize(nbRows, false);
3257 vecF.clear();
3258
3259 // get number of integration points
3260 nbIntegrationPts = getGaussPtsSlave().size2();
3261
3262 // integrate local matrix for entity block
3263 CHKERR iNtegrate(row_data);
3264
3265 // assemble local matrix
3266 CHKERR aSsemble(row_data);
3267
3269}
3270
3272SimpleContactProblem::OpCalMatForcesALESlave::iNtegrate(EntData &data) {
3274
3275 int nb_base_fun_col = nbRows / 3;
3276
3277 auto get_tensor_vec = [](VectorDouble &n, const int r) {
3278 return FTensor::Tensor1<double *, 3>(&n(r + 0), &n(r + 1), &n(r + 2));
3279 };
3280
3283 auto lagrange_slave =
3284 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
3285
3286 auto t_F = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->FMat);
3287 auto normal_at_gp =
3288 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr), 0);
3289 auto t_w = getFTensor0IntegrationWeightSlave();
3290 double &area_slave = commonDataSimpleContact->areaSlave;
3291 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
3292 double val_s = t_w * area_slave;
3293 FTensor::Tensor0<double *> t_base_master(&data.getN()(gg, 0));
3294 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
3295 const double s = val_s * t_base_master * lagrange_slave;
3296 auto t_assemble_s = get_tensor_vec(vecF, 3 * bbc);
3297 t_assemble_s(i) -= s * t_F(j, i) * normal_at_gp(j);
3298 ++t_base_master;
3299 }
3300 ++t_F;
3301 ++lagrange_slave;
3302 ++t_w;
3303 }
3304
3306}
3307
3309SimpleContactProblem::OpCalMatForcesALESlave::aSsemble(EntData &row_data) {
3311
3312 // get pointer to first global index on row
3313 const int *row_indices = &*row_data.getIndices().data().begin();
3314 auto &data = *commonDataSimpleContact;
3315 if (data.forcesOnlyOnEntitiesRow.empty())
3317
3318 rowIndices.resize(nbRows, false);
3319 noalias(rowIndices) = row_data.getIndices();
3320 row_indices = &rowIndices[0];
3321 VectorDofs &dofs = row_data.getFieldDofs();
3322 VectorDofs::iterator dit = dofs.begin();
3323 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
3324 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
3325 data.forcesOnlyOnEntitiesRow.end()) {
3326 rowIndices[ii] = -1;
3327 }
3328 }
3329
3330 CHKERR VecSetValues(getSNESf(), nbRows, row_indices, &*vecF.data().begin(),
3331 ADD_VALUES);
3332
3334}
3335
3338 EntData &data) {
3340
3341 if (data.getFieldData().size() == 0)
3343
3344 if (type != MBVERTEX)
3346
3350
3351 commonDataSimpleContact->normalVectorSlavePtr->resize(3, false);
3352
3353 auto get_tensor_vec = [](VectorDouble &n) {
3354 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3355 };
3356
3357 auto normal_original_slave =
3358 get_tensor_vec(*commonDataSimpleContact->normalVectorSlavePtr);
3359
3360 commonDataSimpleContact->tangentOneVectorSlavePtr->resize(3, false);
3361 commonDataSimpleContact->tangentOneVectorSlavePtr->clear();
3362
3363 commonDataSimpleContact->tangentTwoVectorSlavePtr->resize(3, false);
3364 commonDataSimpleContact->tangentTwoVectorSlavePtr->clear();
3365
3366 auto tangent_0_slave =
3367 get_tensor_vec(*commonDataSimpleContact->tangentOneVectorSlavePtr);
3368 auto tangent_1_slave =
3369 get_tensor_vec(*commonDataSimpleContact->tangentTwoVectorSlavePtr);
3370
3371 auto t_N = data.getFTensor1DiffN<2>(0, 0);
3372 auto t_dof = data.getFTensor1FieldData<3>();
3373
3374 for (unsigned int dd = 0; dd != 3; ++dd) {
3375 tangent_0_slave(i) += t_dof(i) * t_N(0);
3376 tangent_1_slave(i) += t_dof(i) * t_N(1);
3377 ++t_dof;
3378 ++t_N;
3379 }
3380
3381 normal_original_slave(i) =
3382 FTensor::levi_civita(i, j, k) * tangent_0_slave(j) * tangent_1_slave(k);
3383
3384 const double normal_length =
3385 sqrt(normal_original_slave(i) * normal_original_slave(i));
3386 normal_original_slave(i) = normal_original_slave(i) / normal_length;
3387
3388 commonDataSimpleContact->areaSlave = 0.5 * normal_length;
3389
3391}
3392
3395 EntData &data) {
3397
3398 if (data.getFieldData().size() == 0)
3400
3401 if (type != MBVERTEX)
3403
3407
3408 commonDataSimpleContact->normalVectorMasterPtr->resize(3, false);
3409
3410 auto get_tensor_vec = [](VectorDouble &n) {
3411 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3412 };
3413
3414 auto normal_original_master =
3415 get_tensor_vec(*commonDataSimpleContact->normalVectorMasterPtr);
3416
3417 commonDataSimpleContact->tangentOneVectorMasterPtr->resize(3, false);
3418 commonDataSimpleContact->tangentOneVectorMasterPtr->clear();
3419
3420 commonDataSimpleContact->tangentTwoVectorMasterPtr->resize(3, false);
3421 commonDataSimpleContact->tangentTwoVectorMasterPtr->clear();
3422
3423 auto tangent_0_master =
3424 get_tensor_vec(*commonDataSimpleContact->tangentOneVectorMasterPtr);
3425 auto tangent_1_master =
3426 get_tensor_vec(*commonDataSimpleContact->tangentTwoVectorMasterPtr);
3427
3428 auto t_N = data.getFTensor1DiffN<2>(0, 0);
3429 auto t_dof = data.getFTensor1FieldData<3>();
3430
3431 for (unsigned int dd = 0; dd != 3; ++dd) {
3432 tangent_0_master(i) += t_dof(i) * t_N(0);
3433 tangent_1_master(i) += t_dof(i) * t_N(1);
3434 ++t_dof;
3435 ++t_N;
3436 }
3437
3438 normal_original_master(i) =
3439 FTensor::levi_civita(i, j, k) * tangent_0_master(j) * tangent_1_master(k);
3440
3441 const double normal_length =
3442 sqrt(normal_original_master(i) * normal_original_master(i));
3443 normal_original_master(i) = normal_original_master(i) / normal_length;
3444
3445 commonDataSimpleContact->areaMaster = 0.5 * normal_length;
3446
3448}
3449
3451 EntData &row_data, EntData &col_data) {
3453
3454 // Both sides are needed since both sides contribute their shape
3455 // function to the stiffness matrix
3456 const int nb_row = row_data.getIndices().size();
3457 if (!nb_row)
3459 const int nb_col = col_data.getIndices().size();
3460 if (!nb_col)
3462 const int nb_gauss_pts = row_data.getN().size1();
3463
3464 int nb_base_fun_row = row_data.getFieldData().size() / 3;
3465 int nb_base_fun_col = col_data.getFieldData().size() / 3;
3466
3467 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
3469 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
3470 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
3471 &m(r + 2, c + 2));
3472 };
3473
3474 auto get_tensor_vec = [](VectorDouble &n) {
3475 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3476 };
3477
3478 auto get_tensor_vec_3 = [&](VectorDouble3 &n) {
3479 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3480 };
3481
3485
3486 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
3488 t_n(i, j) = 0;
3489 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
3490 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
3491 return t_n;
3492 };
3493
3494 matLhs.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
3495 matLhs.clear();
3496
3497 auto lagrange_slave =
3498 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
3499
3500 auto t_1 = get_tensor_vec(*commonDataSimpleContact->tangentOneVectorSlavePtr);
3501 auto t_2 = get_tensor_vec(*commonDataSimpleContact->tangentTwoVectorSlavePtr);
3502 auto t_w = getFTensor0IntegrationWeightSlave();
3503
3504 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
3505 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
3506
3507 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
3508
3509 FTensor::Tensor0<double *> t_base_slave(&row_data.getN()(gg, 0));
3510
3511 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
3512 const double s = t_w * lagrange_slave * t_base_slave;
3513
3514 auto t_d_n = make_vec_der(t_N, t_1, t_2);
3515
3516 auto t_assemble_s = get_tensor_from_mat(matLhs, 3 * bbr, 3 * bbc);
3517
3518 t_assemble_s(i, j) += 0.5 * s * t_d_n(i, j);
3519
3520 ++t_base_slave; // update rows
3521 }
3522 ++t_N;
3523 }
3524 ++lagrange_slave;
3525 ++t_w;
3526 }
3528}
3529
3531 EntData &row_data, EntData &col_data) {
3533
3534 // Both sides are needed since both sides contribute their shape
3535 // function to the stiffness matrix
3536 const int nb_row = row_data.getIndices().size();
3537 if (!nb_row)
3539 const int nb_col = col_data.getIndices().size();
3540 if (!nb_col)
3542 const int nb_gauss_pts = row_data.getN().size1();
3543
3544 int nb_base_fun_row = row_data.getFieldData().size() / 3;
3545 int nb_base_fun_col = col_data.getFieldData().size() / 3;
3546
3547 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
3549 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
3550 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
3551 &m(r + 2, c + 2));
3552 };
3553
3554 auto get_tensor_vec = [](VectorDouble &n) {
3555 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3556 };
3557
3558 auto get_tensor_vec_3 = [&](VectorDouble3 &n) {
3559 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3560 };
3561
3565
3566 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
3568 t_n(i, j) = 0;
3569 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
3570 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
3571 return t_n;
3572 };
3573
3574 auto lagrange_slave =
3575 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
3576
3577 const double area_s = commonDataSimpleContact->areaSlave;
3578
3579 const double area_m = commonDataSimpleContact->areaMaster;
3580
3581 auto t_1 = get_tensor_vec(*commonDataSimpleContact->tangentOneVectorSlavePtr);
3582 auto t_2 = get_tensor_vec(*commonDataSimpleContact->tangentTwoVectorSlavePtr);
3583
3584 auto t_const_unit_slave =
3585 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
3586
3587 auto t_w = getFTensor0IntegrationWeightMaster();
3588 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
3589
3590 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
3591
3592 const double mult_s = 0.5 * t_w * lagrange_slave * area_m / area_s;
3593
3594 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
3595
3596 FTensor::Tensor0<double *> t_base_master(&row_data.getN()(gg, 0));
3597
3598 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
3599 const double s = mult_s * t_base_master;
3600
3601 auto t_d_n = make_vec_der(t_N, t_1, t_2);
3602
3603 auto t_assemble_s = get_tensor_from_mat(matLhs, 3 * bbr, 3 * bbc);
3604
3605 t_assemble_s(i, j) -=
3606 s * (-t_const_unit_slave(i) * t_d_n(k, j) * t_const_unit_slave(k) +
3607 t_d_n(i, j));
3608
3609 ++t_base_master; // update rows
3610 }
3611 ++t_N;
3612 }
3613 ++lagrange_slave;
3614 ++t_w;
3615 }
3616
3618}
3619
3622 EntData &row_data, EntData &col_data) {
3624
3625 // Both sides are needed since both sides contribute their shape
3626 // function to the stiffness matrix
3627 const int nb_row = row_data.getIndices().size();
3628 if (!nb_row)
3630 const int nb_col = col_data.getIndices().size();
3631 if (!nb_col)
3633 const int nb_gauss_pts = row_data.getN().size1();
3634
3635 int nb_base_fun_row = row_data.getFieldData().size() / 3;
3636 int nb_base_fun_col = col_data.getFieldData().size() / 3;
3637
3638 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
3640 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
3641 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
3642 &m(r + 2, c + 2));
3643 };
3644
3645 auto get_tensor_vec = [](VectorDouble &n) {
3646 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3647 };
3648
3649 auto get_tensor_vec_3 = [&](VectorDouble3 &n) {
3650 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3651 };
3652
3656
3657 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
3659 t_n(i, j) = 0;
3660 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
3661 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
3662 return t_n;
3663 };
3664
3665 matLhs.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
3666 matLhs.clear();
3667
3668 auto lagrange_slave =
3669 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
3670
3671 auto t_1 =
3672 get_tensor_vec(*commonDataSimpleContact->tangentOneVectorMasterPtr);
3673 auto t_2 =
3674 get_tensor_vec(*commonDataSimpleContact->tangentTwoVectorMasterPtr);
3675
3676 auto t_w = getFTensor0IntegrationWeightMaster();
3677
3678 auto t_const_unit_master =
3679 get_tensor_vec(*(commonDataSimpleContact->normalVectorMasterPtr));
3680
3681 auto t_const_unit_slave =
3682 get_tensor_vec(*(commonDataSimpleContact->normalVectorSlavePtr));
3683
3684 const double area_m = commonDataSimpleContact->areaMaster;
3685 const double area_s = commonDataSimpleContact->areaSlave;
3686
3687 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
3688
3689 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
3690
3691 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
3692 const double mult_m = 0.5 * t_w * lagrange_slave;
3693 FTensor::Tensor0<double *> t_base_master(&row_data.getN()(gg, 0));
3694
3695 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
3696 const double s = mult_m * t_base_master;
3697
3698 auto t_d_n = make_vec_der(t_N, t_1, t_2);
3699
3700 auto t_assemble_s = get_tensor_from_mat(matLhs, 3 * bbr, 3 * bbc);
3701
3702 t_assemble_s(i, j) -=
3703 s * t_d_n(k, j) * t_const_unit_master(k) * t_const_unit_slave(i);
3704
3705 ++t_base_master; // update rows
3706 }
3707 ++t_N;
3708 }
3709 ++lagrange_slave;
3710 ++t_w;
3711 }
3712
3714}
3715
3717 EntData &row_data, EntData &col_data) {
3719
3720 const int nb_row = row_data.getIndices().size();
3721 if (!nb_row)
3723 const int nb_col = col_data.getIndices().size();
3724 if (!nb_col)
3726 const int nb_gauss_pts = row_data.getN().size1();
3727 int nb_base_fun_row = row_data.getFieldData().size();
3728 int nb_base_fun_col = col_data.getFieldData().size() / 3;
3729
3730 matLhs.resize(nb_base_fun_row, 3 * nb_base_fun_col, false);
3731 matLhs.clear();
3732
3733 auto get_tensor_vec = [](VectorDouble &n) {
3734 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3735 };
3736
3737 auto get_vec_from_mat = [](MatrixDouble &m, const int r, const int c) {
3738 return FTensor::Tensor1<double *, 3>(&m(r + 0, c + 0), &m(r + 0, c + 1),
3739 &m(r + 0, c + 2));
3740 };
3741
3745
3746 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
3748 t_n(i, j) = 0;
3749 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
3750 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
3751 return t_n;
3752 };
3753
3754 auto x_m = getFTensor1FromMat<3>(
3755 *commonDataSimpleContact->positionAtGaussPtsMasterPtr);
3756 auto x_s = getFTensor1FromMat<3>(
3757 *commonDataSimpleContact->positionAtGaussPtsSlavePtr);
3758 auto t_lagrange_slave =
3759 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
3760
3761 const double length_normal = commonDataSimpleContact->areaSlave;
3762
3763 auto normal_at_gp =
3764 get_tensor_vec(*commonDataSimpleContact->normalVectorSlavePtr);
3765
3766 auto t_w = getFTensor0IntegrationWeightSlave();
3767 auto t_1 = get_tensor_vec(*commonDataSimpleContact->tangentOneVectorSlavePtr);
3768 auto t_2 = get_tensor_vec(*commonDataSimpleContact->tangentTwoVectorSlavePtr);
3769 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
3770 const double cn_value = *cNPtr.get();
3771 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
3772 double val_s = t_w * 0.5;
3773 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
3774
3775 for (int bbc = 0; bbc != nb_base_fun_col; ++bbc) {
3776 FTensor::Tensor0<double *> t_base_lambda(&row_data.getN()(gg, 0));
3777
3778 for (int bbr = 0; bbr != nb_base_fun_row; ++bbr) {
3779 const double s = val_s * t_base_lambda;
3780
3781 auto t_d_n = make_vec_der(t_N, t_1, t_2);
3782
3783 auto assemble_mat = get_vec_from_mat(matLhs, bbr, 3 * bbc);
3784
3785 assemble_mat(j) -=
3786 0.5 *
3787 (t_d_n(i, j) - normal_at_gp(i) * t_d_n(k, j) * normal_at_gp(k)) *
3788 (x_s(i) - x_m(i)) * s * cn_value *
3789 (1 + SimpleContactProblem::Sign(t_lagrange_slave -
3790 cn_value * t_gap_gp));
3791
3792 assemble_mat(j) += t_d_n(i, j) * normal_at_gp(i) * s *
3794 cn_value, t_gap_gp, t_lagrange_slave);
3795
3796 ++t_base_lambda; // update rows
3797 }
3798 ++t_N;
3799 }
3800
3801 ++x_m;
3802 ++x_s;
3803 ++t_lagrange_slave;
3804 ++t_gap_gp;
3805 ++t_w;
3806 }
3807
3809}
3810
3812 int row_side, int col_side, EntityType row_type, EntityType col_type,
3813 EntData &row_data, EntData &col_data) {
3814
3816 if (row_type != MBVERTEX)
3818 row_nb_dofs = row_data.getIndices().size();
3819 if (!row_nb_dofs)
3821 col_nb_dofs = col_data.getIndices().size();
3822 if (!col_nb_dofs)
3824 nb_gauss_pts = row_data.getN().size1();
3825
3826 nb_base_fun_row = row_data.getFieldData().size() / rankRow;
3827 nb_base_fun_col = col_data.getFieldData().size() / rankCol;
3828
3829 matLhs.resize(rankRow * nb_base_fun_row, rankCol * nb_base_fun_col, false);
3830 matLhs.clear();
3831
3832 // integrate local matrix for entity block
3833 CHKERR iNtegrate(row_data, col_data);
3834
3835 // assemble local matrix
3836 CHKERR aSsemble(row_data, col_data);
3837
3839}
3840
3842 int row_side, int col_side, EntityType row_type, EntityType col_type,
3843 EntData &row_data, EntData &col_data) {
3844
3846 row_nb_dofs = row_data.getIndices().size();
3847 if (!row_nb_dofs)
3849 col_nb_dofs = col_data.getIndices().size();
3850 if (!col_nb_dofs)
3852 nb_gauss_pts = row_data.getN().size1();
3853
3854 nb_base_fun_row = row_data.getFieldData().size() / rankRow;
3855 nb_base_fun_col = col_data.getFieldData().size() / rankCol;
3856
3857 matLhs.resize(rankRow * nb_base_fun_row, rankCol * nb_base_fun_col, false);
3858 matLhs.clear();
3859
3860 // integrate local matrix for entity block
3861 CHKERR iNtegrate(row_data, col_data);
3862
3863 // assemble local matrix
3864 CHKERR aSsemble(row_data, col_data);
3865
3867}
3868
3870 EntityType type,
3871 EntData &data) {
3873
3874 if (data.getIndices().size() == 0)
3876
3877 const int nb_gauss_pts = data.getN().size1();
3878
3879 vecR.resize(CommonDataSimpleContact::LAST_ELEMENT, false);
3880 vecR.clear();
3881
3882 commonDataSimpleContact->gaussPtsStatePtr->resize(nb_gauss_pts, false);
3883 commonDataSimpleContact->gaussPtsStatePtr->clear();
3884
3885 auto t_state_gp =
3886 getFTensor0FromVec(*commonDataSimpleContact->gaussPtsStatePtr);
3887
3888 auto t_lagrange_slave =
3889 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
3890 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
3891
3892 for (int gg = 0; gg != nb_gauss_pts; gg++) {
3894
3895 if (!almFlag &&
3896 SimpleContactProblem::State(cN, t_gap_gp, t_lagrange_slave)) {
3898 t_state_gp = 1;
3899 }
3900
3901 if (almFlag &&
3902 SimpleContactProblem::StateALM(cN, t_gap_gp, t_lagrange_slave)) {
3904 t_state_gp = 1;
3905 }
3906
3907 ++t_lagrange_slave;
3908 ++t_gap_gp;
3909 ++t_state_gp;
3910 } // for gauss points
3911
3912 constexpr std::array<int, 2> indices = {CommonDataSimpleContact::ACTIVE,
3914 CHKERR VecSetValues(commonDataSimpleContact->gaussPtsStateVec, 2,
3915 indices.data(), &vecR[0], ADD_VALUES);
3916
3918}
3919
3922 EntData &col_data) {
3923
3925
3926 // get pointer to first global index on row
3927 const int *row_indices = &*row_data.getIndices().data().begin();
3928 // get pointer to first global index on column
3929 const int *col_indices = &*col_data.getIndices().data().begin();
3930
3931 auto &data = *commonDataSimpleContact;
3932 if (data.forcesOnlyOnEntitiesRow.empty())
3934
3935 if (!data.forcesOnlyOnEntitiesRow.empty()) {
3936 rowIndices.resize(row_nb_dofs, false);
3937 noalias(rowIndices) = row_data.getIndices();
3938 row_indices = &rowIndices[0];
3939 VectorDofs &dofs = row_data.getFieldDofs();
3940 VectorDofs::iterator dit = dofs.begin();
3941 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
3942 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
3943 data.forcesOnlyOnEntitiesRow.end()) {
3944 rowIndices[ii] = -1;
3945 }
3946 }
3947 }
3948
3949 // assemble local matrix
3950 CHKERR MatSetValues(getSNESB(), row_nb_dofs, row_indices, col_nb_dofs,
3951 col_indices, &*matLhs.data().begin(), ADD_VALUES);
3952
3954}
3955
3958 EntData &col_data) {
3959
3961
3962 // get pointer to first global index on row
3963 const int *row_indices = &*row_data.getIndices().data().begin();
3964 // get pointer to first global index on column
3965 const int *col_indices = &*col_data.getIndices().data().begin();
3966
3967 // assemble local matrix
3968 CHKERR MatSetValues(getSNESB(), row_nb_dofs, row_indices, col_nb_dofs,
3969 col_indices, &*matLhs.data().begin(), ADD_VALUES);
3971}
3972
3975 EntData &row_data, EntData &col_data) {
3976
3978
3982
3983 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
3985 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
3986 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
3987 &m(r + 2, c + 2));
3988 };
3989
3990 auto get_tensor_vec = [](VectorDouble &n) {
3991 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
3992 };
3993
3994 auto t_w = getFTensor0IntegrationWeight();
3995
3996 auto t_inv_H = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->invHMat);
3997
3998 auto lagrange_slave =
3999 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
4000
4001 auto normal_at_gp = get_tensor_vec(*normalVector);
4002
4003 const double area_m = aRea;
4004
4005 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
4006
4007 double a = -t_w * lagrange_slave * area_m;
4008
4009 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
4010
4011 int bbc = 0;
4012 for (; bbc != nb_base_fun_col; ++bbc) {
4013
4014 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
4015
4016 int bbr = 0;
4017 for (; bbr != nb_base_fun_row; ++bbr) {
4018
4019 auto t_assemble = get_tensor2(matLhs, 3 * bbr, 3 * bbc);
4020 // TODO: handle hoGeometry
4021
4022 t_assemble(i, j) += a * t_row_base * t_inv_H(k, i) *
4023 t_col_diff_base(k) * normal_at_gp(j);
4024
4025 ++t_row_base;
4026 }
4027 ++t_col_diff_base;
4028 }
4029 ++t_w;
4030 ++t_inv_H;
4031 ++lagrange_slave;
4032 }
4033
4035}
4036
4039 EntData &row_data, EntData &col_data) {
4040
4042
4048
4049 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
4051 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
4052 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
4053 &m(r + 2, c + 2));
4054 };
4055
4056 auto get_tensor_vec = [](VectorDouble &n) {
4057 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
4058 };
4059
4060 auto t_w = getFTensor0IntegrationWeight();
4061
4062 auto t_h = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->hMat);
4063 auto t_inv_H = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->invHMat);
4064
4065 auto lagrange_slave =
4066 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
4067
4068 auto normal_at_gp = get_tensor_vec(*normalVector);
4069
4070 const double area_m = aRea;
4071
4072 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
4073
4074 const double a = -t_w * lagrange_slave * area_m;
4075 auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
4076
4077 int bbc = 0;
4078 for (; bbc != nb_base_fun_col; ++bbc) {
4079
4080 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
4081
4082 int bbr = 0;
4083 for (; bbr != nb_base_fun_row; ++bbr) {
4084
4085 auto t_assemble = get_tensor2(matLhs, 3 * bbr, 3 * bbc);
4086
4087 t_assemble(i, j) -= a * t_row_base * t_inv_H(l, j) *
4088 t_col_diff_base(m) * t_inv_H(m, i) * t_h(k, l) *
4089 normal_at_gp(k);
4090
4091 ++t_row_base;
4092 }
4093 ++t_col_diff_base;
4094 }
4095 ++t_w;
4096 ++t_h;
4097 ++t_inv_H;
4098 ++lagrange_slave;
4099 }
4100
4102}
4103
4106 EntData &row_data, EntData &col_data) {
4107
4109
4113
4114 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
4116 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
4117 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
4118 &m(r + 2, c + 2));
4119 };
4120
4121 auto get_tensor_vec = [](VectorDouble &n) {
4122 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
4123 };
4124
4125 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
4127 t_n(i, j) = 0;
4128 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
4129 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
4130 return t_n;
4131 };
4132
4133 auto t_F = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->FMat);
4134
4135 auto t_w = getFTensor0IntegrationWeightMaster();
4136 auto lagrange_slave =
4137 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
4138 auto t_1 =
4139 get_tensor_vec(*commonDataSimpleContact->tangentOneVectorMasterPtr);
4140 auto t_2 =
4141 get_tensor_vec(*commonDataSimpleContact->tangentTwoVectorMasterPtr);
4142 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
4143
4144 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
4145 const double val = 0.5 * t_w * lagrange_slave;
4146
4147 int bbc = 0;
4148 for (; bbc != nb_base_fun_col; ++bbc) {
4149
4150 FTensor::Tensor0<double *> t_base(&row_data.getN()(gg, 0));
4151
4152 int bbr = 0;
4153 for (; bbr != nb_base_fun_row; ++bbr) {
4154
4155 auto t_assemble = get_tensor2(matLhs, 3 * bbr, 3 * bbc);
4156 // TODO: handle hoGeometry
4157 auto t_d_n = make_vec_der(t_N, t_1, t_2);
4158 t_assemble(i, k) -= val * t_base * t_F(j, i) * t_d_n(j, k);
4159
4160 ++t_base;
4161 }
4162 ++t_N;
4163 }
4164 ++t_F;
4165 ++t_w;
4166 ++lagrange_slave;
4167 }
4168
4170}
4171
4174 EntData &row_data, EntData &col_data) {
4175
4177
4181
4182 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
4184 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
4185 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
4186 &m(r + 2, c + 2));
4187 };
4188
4189 auto get_tensor_vec = [](VectorDouble &n) {
4190 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
4191 };
4192
4193 auto make_vec_der = [&](auto t_N, auto t_1, auto t_2) {
4195 t_n(i, j) = 0;
4196 t_n(i, j) += FTensor::levi_civita(i, j, k) * t_2(k) * t_N(0);
4197 t_n(i, j) -= FTensor::levi_civita(i, j, k) * t_1(k) * t_N(1);
4198 return t_n;
4199 };
4200
4201 // commonDataSimpleContact->faceRowData = nullptr;
4202
4203 auto t_F = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->FMat);
4204
4205 auto t_w = getFTensor0IntegrationWeightMaster();
4206 auto lagrange_slave =
4207 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
4208 auto t_1 = get_tensor_vec(*commonDataSimpleContact->tangentOneVectorSlavePtr);
4209 auto t_2 = get_tensor_vec(*commonDataSimpleContact->tangentTwoVectorSlavePtr);
4210 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
4211
4212 auto t_N = col_data.getFTensor1DiffN<2>(gg, 0);
4213 const double val = 0.5 * t_w * lagrange_slave;
4214
4215 int bbc = 0;
4216 for (; bbc != nb_base_fun_col; ++bbc) {
4217
4218 FTensor::Tensor0<double *> t_base(&row_data.getN()(gg, 0));
4219
4220 int bbr = 0;
4221 for (; bbr != nb_base_fun_row; ++bbr) {
4222
4223 auto t_assemble = get_tensor2(matLhs, 3 * bbr, 3 * bbc);
4224 // TODO: handle hoGeometry
4225 auto t_d_n = make_vec_der(t_N, t_1, t_2);
4226 t_assemble(i, k) -= val * t_base * t_F(j, i) * t_d_n(j, k);
4227
4228 ++t_base;
4229 }
4230 ++t_N;
4231 }
4232 ++t_F;
4233 ++t_w;
4234 ++lagrange_slave;
4235 }
4236
4238}
4239
4243 EntitiesFieldData::EntData &col_data) {
4244
4246
4249
4250 auto get_tensor_vec = [](VectorDouble &n) {
4251 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
4252 };
4253
4254 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
4255 return FTensor::Tensor1<double *, 3>(&m(r + 0, c + 0), &m(r + 1, c + 0),
4256 &m(r + 2, c + 0));
4257 };
4258
4259 auto t_F = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->FMat);
4260
4261 auto t_w = getFTensor0IntegrationWeightMaster();
4262
4263 auto normal_master_at_gp =
4264 get_tensor_vec(*commonDataSimpleContact->normalVectorMasterPtr);
4265
4266 const double area_m = commonDataSimpleContact->areaMaster;
4267
4268 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
4269
4270 FTensor::Tensor0<double *> t_col_base(&col_data.getN()(gg, 0));
4271
4272 const double val = t_w * area_m;
4273
4274 int bbc = 0;
4275 for (; bbc != nb_base_fun_col; ++bbc) {
4276
4277 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
4278
4279 int bbr = 0;
4280 for (; bbr != nb_base_fun_row; ++bbr) {
4281
4282 auto t_assemble = get_tensor_from_mat(matLhs, 3 * bbr, bbc);
4283
4284 t_assemble(i) -=
4285 val * t_row_base * t_F(j, i) * normal_master_at_gp(j) * t_col_base;
4286
4287 ++t_row_base;
4288 }
4289 ++t_col_base;
4290 }
4291 ++t_F;
4292 ++t_w;
4293 }
4294
4296}
4297
4301 EntitiesFieldData::EntData &col_data) {
4302
4304
4307
4308 auto get_tensor_vec = [](VectorDouble &n) {
4309 return FTensor::Tensor1<double *, 3>(&n(0), &n(1), &n(2));
4310 };
4311
4312 auto get_tensor_from_mat = [](MatrixDouble &m, const int r, const int c) {
4313 return FTensor::Tensor1<double *, 3>(&m(r + 0, c + 0), &m(r + 1, c + 0),
4314 &m(r + 2, c + 0));
4315 };
4316
4317 auto t_F = getFTensor2FromMat<3, 3>(*commonDataSimpleContact->FMat);
4318
4319 auto t_w = getFTensor0IntegrationWeightSlave();
4320
4321 auto normal_master_at_gp =
4322 get_tensor_vec(*commonDataSimpleContact->normalVectorSlavePtr);
4323
4324 const double area_m = commonDataSimpleContact->areaSlave;
4325
4326 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
4327
4328 FTensor::Tensor0<double *> t_col_base(&col_data.getN()(gg, 0));
4329
4330 const double val = t_w * area_m;
4331
4332 int bbc = 0;
4333 for (; bbc != nb_base_fun_col; ++bbc) {
4334
4335 FTensor::Tensor0<double *> t_row_base(&row_data.getN()(gg, 0));
4336
4337 int bbr = 0;
4338 for (; bbr != nb_base_fun_row; ++bbr) {
4339
4340 auto t_assemble = get_tensor_from_mat(matLhs, 3 * bbr, bbc);
4341
4342 t_assemble(i) -=
4343 val * t_row_base * t_F(j, i) * normal_master_at_gp(j) * t_col_base;
4344
4345 ++t_row_base;
4346 }
4347 ++t_col_base;
4348 }
4349 ++t_F;
4350 ++t_w;
4351 }
4352
4354}
4355
4357 int row_side, int col_side, EntityType row_type, EntityType col_type,
4358 EntData &row_data, EntData &col_data) {
4359
4361 if (row_type != MBVERTEX)
4363 row_nb_dofs = row_data.getIndices().size();
4364 if (!row_nb_dofs)
4366 col_nb_dofs = col_data.getIndices().size();
4367 if (!col_nb_dofs)
4369 nb_gauss_pts = row_data.getN().size1();
4370
4371 nb_base_fun_row = row_data.getFieldData().size() / 3;
4372 nb_base_fun_col = col_data.getFieldData().size() / 3;
4373
4374 matLhs.resize(3 * nb_base_fun_row, 3 * nb_base_fun_col, false);
4375 matLhs.clear();
4376
4377 normalVector->resize(3, false);
4378 tangentOne->resize(3, false);
4379 tangentTwo->resize(3, false);
4380
4381 if (isMaster) {
4382 normalVector = commonDataSimpleContact->normalVectorMasterPtr;
4383 tangentOne = commonDataSimpleContact->tangentOneVectorMasterPtr;
4384 tangentTwo = commonDataSimpleContact->tangentOneVectorMasterPtr;
4385 aRea = commonDataSimpleContact->areaMaster;
4386 } else {
4387 normalVector = commonDataSimpleContact->normalVectorSlavePtr;
4388 tangentOne = commonDataSimpleContact->tangentOneVectorSlavePtr;
4389 tangentTwo = commonDataSimpleContact->tangentOneVectorSlavePtr;
4390 aRea = commonDataSimpleContact->areaSlave;
4391 }
4392
4393 // integrate local matrix for entity block
4394 CHKERR iNtegrate(row_data, col_data);
4395
4396 // assemble local matrix
4397 CHKERR aSsemble(row_data, col_data);
4398
4400}
4401
4403 EntData &row_data, EntData &col_data) {
4404
4406
4407 // get pointer to first global index on row
4408 const int *row_indices = &*row_data.getIndices().data().begin();
4409 // get pointer to first global index on column
4410 const int *col_indices = &*col_data.getIndices().data().begin();
4411
4412 auto &data = *commonDataSimpleContact;
4413 if (!data.forcesOnlyOnEntitiesRow.empty()) {
4414 rowIndices.resize(row_nb_dofs, false);
4415 noalias(rowIndices) = row_data.getIndices();
4416 row_indices = &rowIndices[0];
4417 VectorDofs &dofs = row_data.getFieldDofs();
4418 VectorDofs::iterator dit = dofs.begin();
4419 for (int ii = 0; dit != dofs.end(); ++dit, ++ii) {
4420 if (data.forcesOnlyOnEntitiesRow.find((*dit)->getEnt()) ==
4421 data.forcesOnlyOnEntitiesRow.end()) {
4422 rowIndices[ii] = -1;
4423 }
4424 }
4425 }
4426
4427 // assemble local matrix
4428 CHKERR MatSetValues(getSNESB(), row_nb_dofs, row_indices, col_nb_dofs,
4429 col_indices, &*matLhs.data().begin(), ADD_VALUES);
4430
4432}
4433
4434// setup operators for calculation of active set
4436 boost::shared_ptr<SimpleContactElement> fe_rhs_simple_contact_ale,
4437 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
4438 const string field_name, const string mesh_node_field_name,
4439 const string lagrange_field_name, const string side_fe_name) {
4441
4442 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
4443 fe_mat_side_rhs_master = boost::make_shared<
4445
4446 fe_mat_side_rhs_master->getOpPtrVector().push_back(
4448 mesh_node_field_name, common_data_simple_contact->HMat));
4449 fe_mat_side_rhs_master->getOpPtrVector().push_back(
4451 field_name, common_data_simple_contact->hMat));
4452
4453 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
4454 fe_mat_side_rhs_slave = boost::make_shared<
4456
4457 fe_mat_side_rhs_slave->getOpPtrVector().push_back(
4459 mesh_node_field_name, common_data_simple_contact->HMat));
4460 fe_mat_side_rhs_slave->getOpPtrVector().push_back(
4462 field_name, common_data_simple_contact->hMat));
4463
4464 fe_rhs_simple_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
4465 "MESH_NODE_POSITIONS", common_data_simple_contact));
4466
4467 fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4468 new OpGetNormalMasterALE("MESH_NODE_POSITIONS",
4469 common_data_simple_contact));
4470
4471 // fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4472 // new OpGetPositionAtGaussPtsMaster(field_name,
4473 // common_data_simple_contact));
4474
4475 // fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4476 // new OpGetPositionAtGaussPtsSlave(field_name,
4477 // common_data_simple_contact));
4478
4479 // fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4480 // new OpGetGapSlave(field_name, common_data_simple_contact));
4481
4482 fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4483 new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
4484 common_data_simple_contact));
4485
4486 fe_mat_side_rhs_master->getOpPtrVector().push_back(new OpCalculateDeformation(
4487 mesh_node_field_name, common_data_simple_contact, false));
4488
4489 fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4490 new OpLoopForSideOfContactPrism(mesh_node_field_name,
4491 fe_mat_side_rhs_master, side_fe_name,
4492 ContactOp::FACEMASTER));
4493
4494 fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4495 new OpCalMatForcesALEMaster(mesh_node_field_name,
4496 common_data_simple_contact));
4497
4498 fe_mat_side_rhs_slave->getOpPtrVector().push_back(new OpCalculateDeformation(
4499 mesh_node_field_name, common_data_simple_contact, false));
4500
4501 fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4502 new OpLoopForSideOfContactPrism(mesh_node_field_name,
4503 fe_mat_side_rhs_slave, side_fe_name,
4504 ContactOp::FACESLAVE));
4505
4506 fe_rhs_simple_contact_ale->getOpPtrVector().push_back(
4507 new OpCalMatForcesALESlave(mesh_node_field_name,
4508 common_data_simple_contact));
4510}
4511
4513 boost::shared_ptr<SimpleContactElement> fe_lhs_simple_contact_ale,
4514 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
4515 const string field_name, const string mesh_node_field_name,
4516 const string lagrange_field_name, const string side_fe_name) {
4518
4519 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
4520 mesh_node_field_name, common_data_simple_contact));
4521
4522 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4523 new OpGetNormalMasterALE(mesh_node_field_name,
4524 common_data_simple_contact));
4525
4526 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4527 new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
4528 common_data_simple_contact));
4529
4530 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
4531 feMatSideLhs_dx = boost::make_shared<
4533
4534 feMatSideLhs_dx->getOpPtrVector().push_back(
4536 mesh_node_field_name, common_data_simple_contact->HMat));
4537
4538 feMatSideLhs_dx->getOpPtrVector().push_back(
4540 field_name, common_data_simple_contact->hMat));
4541
4542 // // Master derivative over spatial
4543 feMatSideLhs_dx->getOpPtrVector().push_back(new OpCalculateDeformation(
4544 mesh_node_field_name, common_data_simple_contact, false));
4545
4546 feMatSideLhs_dx->getOpPtrVector().push_back(
4548 mesh_node_field_name, mesh_node_field_name,
4549 common_data_simple_contact, true));
4550
4551 feMatSideLhs_dx->getOpPtrVector().push_back(
4553 mesh_node_field_name, field_name, common_data_simple_contact, true));
4554
4555 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4556 new OpLoopForSideOfContactPrism(mesh_node_field_name, feMatSideLhs_dx,
4557 side_fe_name, ContactOp::FACEMASTER));
4558
4559 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4561 mesh_node_field_name, mesh_node_field_name,
4562 common_data_simple_contact, POSITION_RANK, POSITION_RANK));
4563
4564 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4566 mesh_node_field_name, lagrange_field_name, common_data_simple_contact,
4568
4569 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnContactPrismSide>
4570 feMatSideLhsSlave_dx = boost::make_shared<
4572
4573 feMatSideLhsSlave_dx->getOpPtrVector().push_back(
4575 mesh_node_field_name, common_data_simple_contact->HMat));
4576
4577 feMatSideLhsSlave_dx->getOpPtrVector().push_back(
4579 field_name, common_data_simple_contact->hMat));
4580
4581 feMatSideLhsSlave_dx->getOpPtrVector().push_back(new OpCalculateDeformation(
4582 mesh_node_field_name, common_data_simple_contact, false));
4583
4584 feMatSideLhsSlave_dx->getOpPtrVector().push_back(
4586 mesh_node_field_name, mesh_node_field_name,
4587 common_data_simple_contact, false));
4588
4589 feMatSideLhsSlave_dx->getOpPtrVector().push_back(
4591 mesh_node_field_name, field_name, common_data_simple_contact, false));
4592
4593 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4594 new OpLoopForSideOfContactPrism(mesh_node_field_name,
4595 feMatSideLhsSlave_dx, side_fe_name,
4596 ContactOp::FACESLAVE));
4597
4598 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4600 mesh_node_field_name, mesh_node_field_name,
4601 common_data_simple_contact, POSITION_RANK, POSITION_RANK));
4602
4603 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4605 mesh_node_field_name, lagrange_field_name, common_data_simple_contact,
4607
4609}
4610
4612 boost::shared_ptr<SimpleContactElement> fe_lhs_simple_contact_ale,
4613 boost::shared_ptr<CommonDataSimpleContact> common_data_simple_contact,
4614 const string field_name, const string mesh_node_field_name,
4615 const string lagrange_field_name, bool is_eigen_pos_field,
4616 string eigen_pos_field_name) {
4618
4619 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(new OpGetNormalSlaveALE(
4620 mesh_node_field_name, common_data_simple_contact));
4621
4622 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4623 new OpGetNormalMasterALE(mesh_node_field_name,
4624 common_data_simple_contact));
4625
4626 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4628 common_data_simple_contact));
4629
4630 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4631 new OpGetPositionAtGaussPtsSlave(field_name, common_data_simple_contact));
4632
4633 if (is_eigen_pos_field) {
4634 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4636 eigen_pos_field_name, common_data_simple_contact));
4637
4638 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4640 eigen_pos_field_name, common_data_simple_contact));
4641 }
4642
4643 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4644 new OpGetGapSlave(field_name, common_data_simple_contact));
4645
4646 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4647 new OpGetLagMulAtGaussPtsSlave(lagrange_field_name,
4648 common_data_simple_contact));
4649
4650 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4651 new OpContactTractionSlaveSlave_dX(field_name, mesh_node_field_name,
4652 common_data_simple_contact,
4654
4655 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4656 new OpContactTractionMasterSlave_dX(field_name, mesh_node_field_name,
4657 common_data_simple_contact,
4659
4660 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4661 new OpContactTractionMasterMaster_dX(field_name, mesh_node_field_name,
4662 common_data_simple_contact,
4664
4665 fe_lhs_simple_contact_ale->getOpPtrVector().push_back(
4667 lagrange_field_name, mesh_node_field_name, cnValuePtr,
4668 common_data_simple_contact, LAGRANGE_RANK, POSITION_RANK));
4670}
4671
4673 EntityType type,
4674 EntData &data) {
4676
4677 if (data.getIndices().size() == 0)
4679
4680 if (!postProcSurface.empty()) {
4681 const EntityHandle prism_ent = getFEEntityHandle();
4682 Range tri_ents;
4683 auto &m_field = this->getPtrFE()->mField;
4684 CHKERR m_field.get_moab().get_adjacencies(&prism_ent, 1, 2, false, tri_ents,
4685 moab::Interface::UNION);
4686 tri_ents = tri_ents.subset_by_type(MBTRI);
4687 if (intersect(postProcSurface, tri_ents).empty())
4689 }
4690
4691 const int nb_gauss_pts = data.getN().size1();
4692
4693 int nb_base_fun_col = data.getFieldData().size();
4694 const double area_s =
4695 commonDataSimpleContact->areaSlave; // same area in master and slave
4696
4697 vecR.resize(CommonDataSimpleContact::LAST_ELEMENT, false);
4698 vecR.clear();
4699
4700 auto t_lagrange_slave =
4701 getFTensor0FromVec(*commonDataSimpleContact->lagMultAtGaussPtsPtr);
4702 auto t_gap_gp = getFTensor0FromVec(*commonDataSimpleContact->gapPtr);
4703 auto t_w = getFTensor0IntegrationWeightSlave();
4704
4705 for (int gg = 0; gg != nb_gauss_pts; gg++) {
4706 const double val_s = t_w * area_s;
4707 vecR[CommonDataSimpleContact::TOTAL] += val_s;
4708 bool gap_below_tolerance =
4709 postProcGapTol > std::numeric_limits<double>::epsilon() &&
4710 t_gap_gp < postProcGapTol;
4711 if (gap_below_tolerance) {
4712 vecR[CommonDataSimpleContact::ACTIVE] += val_s;
4713 } else {
4714 if (!almFlag &&
4715 SimpleContactProblem::State(cN, t_gap_gp, t_lagrange_slave)) {
4716 vecR[CommonDataSimpleContact::ACTIVE] += val_s;
4717 }
4718 if (almFlag &&
4719 SimpleContactProblem::StateALM(cN, t_gap_gp, t_lagrange_slave)) {
4720 vecR[CommonDataSimpleContact::ACTIVE] += val_s;
4721 }
4722 }
4723
4724 ++t_lagrange_slave;
4725 ++t_gap_gp;
4726 ++t_w;
4727 } // for gauss points
4728
4729 constexpr std::array<int, 2> indices = {
4732 };
4733 CHKERR VecSetValues(commonDataSimpleContact->contactAreaVec, 2,
4734 indices.data(), &vecR[0], ADD_VALUES);
4735
4737}
static Index< 'J', 3 > J
constexpr double a
static const double eps
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
@ F
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double tol
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
Definition: Templates.hpp:1620
static auto determinantTensor2by2(T &t)
Calculate the determinant of a 2x2 matrix or a tensor of rank 2.
Definition: Templates.hpp:1524
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1559
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1511
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
ublas::vector< FEDofEntity *, DofsAllocator > VectorDofs
constexpr IntegrationType I
constexpr AssemblyType A
constexpr double t
plate stiffness
Definition: plate.cpp:59
constexpr auto field_name
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1FieldData()
Return FTensor of rank 1, i.e. vector from filed data coeffinects.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FieldData()
Resturn scalar files as a FTensor of rank 0.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
static double shapeFunMBTRI1(const double x, const double y)
Definition: Tools.hpp:640
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
static double shapeFunMBTRI0(const double x, const double y)
Definition: Tools.hpp:636
static double shapeFunMBTRI2(const double x, const double y)
Definition: Tools.hpp:644
static MoFEMErrorCode shapeFunMBTRI(double *shape, const double *ksi, const double *eta, const int nb)
Calculate shape functions on triangle.
Definition: Tools.hpp:649
Base volume element used to integrate on contact surface (could be extended to other volume elements)...
RHS-operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrates Lagrange multipliers virtual work on slave surface and assembles components to the RHS glo...
RHS-operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Integrates Lagrange multipliers virtual work on slave surface and assembles components to the RHS glo...
LHS-operator for the simple contact element with Augmented Lagrangian Method.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates virtual work on master side , , derivative with respect to Lagrange multipliers on slave s...
LHS-operator for the simple contact element with Augmented Lagrangian Method.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates virtual work on slave side , , derivative with respect to Lagrange multipliers on slave si...
LHS-operator for the simple contact element with Augmented Lagrangian Method.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrates virtual work on master side , , derivative with respect to spatial positions of the master...
LHS-operator for the simple contact element with Augmented Lagrangian Method.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrates virtual work on master side , derivative with respect to spatial positions on slave side ...
LHS-operator for the simple contact element with Augmented Lagrangian Method.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrates virtual work on slave side, , derivative with respect to spatial positions on master side ...
LHS-operator for the simple contact element with Augmented Lagrangian Method.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrates virtual work on slave side, , derivative with respect to slave spatial positions and assem...
RHS-operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Integrates Lagrange multipliers virtual work on master surface and assembles to global RHS vector.
RHS-operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Integrates Lagrange multipliers virtual work on slave surface and assembles components to the RHS glo...
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates Lagrange multipliers virtual work, derivative with respect to Lagrange multipliers with r...
LHS-operator for the simple contact element.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates and assembles Lagrange multipliers virtual work, derivative with respect to Lagrange mult...
LHS-operator for the simple contact element.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates the complementarity function at slave face gauss points and assembles components to LHS gl...
LHS-operator for the simple contact element.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates linearisation of the complementarity function at slave face gauss points and assembles com...
LHS-operator for the simple contact element.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates linearisation of the complementarity function at slave face gauss points and assembles com...
LHS-operator for the simple contact element.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates linearisation of the complementarity function at slave face gauss points and assembles com...
RHS-operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Integrates the complementarity function at slave face gauss points and assembles components to the RH...
RHS - operator for the contact element (material configuration) Integrates virtual work of contact tr...
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Integrates virtual work of contact traction in the material configuration.
MoFEMErrorCode iNtegrate(EntData &row_data)
MoFEMErrorCode aSsemble(EntData &row_data)
RHS - operator for the contact element (material configuration) Integrates virtual work of contact tr...
Operator for computing deformation gradients in side volumes.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Evaluate gradient of Lagrange multipliers on reference slave surface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Evaluate gradient position on reference master surface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
LHS-operator for the contact element (material configuration)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the contact element (material configuration)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the contact element (material configuration)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the contact element (material configuration)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Compute part of the left-hand side.
LHS-operator for the contact element (material configuration)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Compute part of the left-hand side.
LHS-operator (material configuration) on the side volume of either master or slave side.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates over a face contribution from a side volume.
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
LHS-operator for the simple contact element.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates Lagrange multipliers virtual work on master side, derivative with respect to material pos...
LHS-operator for the simple contact element.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates Lagrange multipliers virtual work on master side, derivative with respect to material pos...
LHS-operator for the simple contact element.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrates Lagrange multipliers virtual work, derivative with respect to material positions on slave...
LHS-operator for the simple contact element.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates the conditions that fulfil KKT conditions at slave face gauss points and assembles compone...
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates the conditions that fulfil KKT conditions at master face gauss points and assembles compon...
LHS-operator for the simple contact element.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Integrates the conditions that fulfil KKT conditions at slave face gauss points and assembles compone...
RHS-operator for the simple contact element for Augmented Lagrangian Method.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for the simple contact element for Augmented Lagrangian Method.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Evaluates gap function at slave face gauss points.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Evaluates unit normal vector to the master surface vector based on material base coordinates.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Evaluates unit normal vector to the master surface vector based on reference base coordinates.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Evaluates unit normal vector to the slave surface vector based on material base coordinates.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Evaluates unit normal vector to the slave surface vector based on reference base coordinates.
Operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Tangent opeerator for contrains equation for change of spatial positions on master and slave.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Calculate tangent operator for contact force for change of integration point positions,...
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Trigers operators for side volume element on slave side for evaluation of the RHS contact traction in...
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for the simple contact element.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
virtual MoFEMErrorCode setGaussPts(int order)
MoFEMErrorCode setContactOperatorsRhs(boost::shared_ptr< SimpleContactElement > fe_rhs_simple_contact, boost::shared_ptr< CommonDataSimpleContact > common_data_simple_contact, string field_name, string lagrange_field_name, bool is_alm=false, bool is_eigen_pos_field=false, string eigen_pos_field_name="EIGEN_SPATIAL_POSITIONS", bool use_reference_coordinates=false)
Function for the simple contact element for C function or ALM approach that sets the user data RHS-op...
MoFEM::Interface & mField
static double ConstrainFunction_dg(const double cn, const double g, const double l)
static constexpr int POSITION_RANK
MoFEMErrorCode setContactOperatorsLhs(boost::shared_ptr< SimpleContactElement > fe_lhs_simple_contact, boost::shared_ptr< CommonDataSimpleContact > common_data_simple_contact, string field_name, string lagrange_field_name, bool is_alm=false, bool is_eigen_pos_field=false, string eigen_pos_field_name="EIGEN_SPATIAL_POSITIONS", bool use_reference_coordinates=false)
Function for the simple contact element for C function or ALM approach that sets the user data LHS-op...
MoFEMErrorCode setMasterForceOperatorsRhs(boost::shared_ptr< SimpleContactElement > fe_lhs_simple_contact, boost::shared_ptr< CommonDataSimpleContact > common_data_simple_contact, string field_name, string lagrange_field_name, bool is_alm=false, bool is_eigen_pos_field=false, string eigen_pos_field_name="EIGEN_SPATIAL_POSITIONS", bool use_reference_coordinates=false)
static double Sign(double x)
Definition:</