56 {
58
59 auto get_material_dofs_from_coords = [&]() {
62 int num_nodes;
65 conn, num_nodes, true);
70 };
71
72 auto get_dofs_data_for_slave_and_master = [&] {
77 for (
size_t n = 0;
n != 3; ++
n) {
78 for (
size_t d = 0;
d != 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(
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);
97 &
fePtr->gaussPtsSlave(1, 0), nb_gauss_pts);
99 &
fePtr->gaussPtsMaster(1, 0), nb_gauss_pts);
101 };
102
104 if (CONVECT_MASTER)
106 else
108 };
109
111 if (CONVECT_MASTER)
113 else
115 };
116
117 auto get_slave_material_coords = [&]() ->
MatrixDouble & {
118 if (CONVECT_MASTER)
120 else
122 };
123
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)
134 else
136 };
137
138 auto get_master_spatial_coords = [&]() ->
MatrixDouble & {
139 if (CONVECT_MASTER)
141 else
143 };
144
146 if (CONVECT_MASTER)
148 else
150 };
151
153 if (CONVECT_MASTER)
155 else
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) {
169 ©_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) {
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 = [&]() {
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 = [&]() {
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) =
277 get_master_n()(gg, 1) =
279 get_master_n()(gg, 2) =
281 };
282
285 double det_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;
303
304 do {
305
306 get_values();
307 assemble();
308 linear_solver();
309 update();
310
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();
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}
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
VectorBoundedArray< double, 3 > VectorDouble3
MoFEMErrorCode determinantTensor2by2(T1 &t, T2 &det)
Calculate determinant 2 by 2.
MoFEMErrorCode invertTensor2by2(T1 &t, T2 &det, T3 &inv_t)
Calculate matrix inverse 2 by 2.
constexpr IntegrationType I
virtual moab::Interface & get_moab()=0
