32 template <
typename T>
inline static T Heav(T a)
33 {
return (a < T(0)) ? T(0) : T(1); }
42 void contact_nonlinear_term::adjust_tensor_size(
void) {
43 sizes_.resize(1); sizes_[0] = 1;
46 case RHS_U_V1:
case RHS_U_V2:
case RHS_U_V4:
47 case RHS_U_V5:
case RHS_U_FRICT_V6:
case RHS_U_FRICT_V7:
48 case RHS_U_FRICT_V8:
case RHS_U_FRICT_V1:
49 case RHS_U_FRICT_V4:
case RHS_U_FRICT_V5:
50 case RHS_L_FRICT_V1:
case RHS_L_FRICT_V2:
case RHS_L_FRICT_V4:
51 case K_UL_V1:
case K_UL_V2:
case K_UL_V3:
52 case UZAWA_PROJ_FRICT:
case UZAWA_PROJ_FRICT_SAXCE:
55 case K_UU_V1:
case K_UU_V2:
56 case K_UL_FRICT_V1:
case K_UL_FRICT_V2:
case K_UL_FRICT_V3:
57 case K_UL_FRICT_V4:
case K_UL_FRICT_V5:
58 case K_UL_FRICT_V7:
case K_UL_FRICT_V8:
59 case K_LL_FRICT_V1:
case K_LL_FRICT_V2:
case K_LL_FRICT_V4:
60 case K_UU_FRICT_V1:
case K_UU_FRICT_V2:
61 case K_UU_FRICT_V3:
case K_UU_FRICT_V4:
case K_UU_FRICT_V5:
62 sizes_.resize(2); sizes_[0] = sizes_[1] = N;
break;
66 lnt.resize(N); lt.resize(N); zt.resize(N); no.resize(N);
67 aux1.resize(1); auxN.resize(N); V.resize(N);
71 void contact_nonlinear_term::friction_law
72 (scalar_type p, scalar_type &tau) {
73 tau = (p > scalar_type(0)) ? tau_adh + f_coeff * p : scalar_type(0);
74 if (tau > tresca_lim) tau = tresca_lim;
77 void contact_nonlinear_term::friction_law
78 (scalar_type p, scalar_type &tau, scalar_type &tau_grad) {
79 if (p <= scalar_type(0)) {
81 tau_grad = scalar_type(0);
84 tau = tau_adh + f_coeff * p;
85 if (tau > tresca_lim) {
87 tau_grad = scalar_type(0);
94 void contact_nonlinear_term::compute
95 (fem_interpolation_context &, bgeot::base_tensor &t) {
97 t.adjust_sizes(sizes_);
98 scalar_type e, augm_ln, rho, rho_grad;
107 t[0] = (ln+gmm::neg(ln-r*(un - g)))/r;
break;
109 t[0] = (un-g) + gmm::pos(ln)/r;
break;
112 t[0] = (Heav(r*(un-g)-ln) - scalar_type(1))/r;
break;
114 t[0] = -Heav(ln)/r;
break;
117 t[0] = -gmm::neg(ln - r*(un - g));
break;
122 t[0] = Heav(un-g - ln);
break;
123 case CONTACT_PRESSURE:
129 for (i=0; i<N; ++i) t[i] = ln * no[i];
132 e = -gmm::neg(ln-r*(un - g));
133 for (i=0; i<N; ++i) t[i] = e * no[i];
137 for (i=0; i<N; ++i) t[i] = e * no[i];
140 e = - gmm::pos(un-g) * r;
141 for (i=0; i<N; ++i) t[i] = e * no[i];
144 e = gmm::neg(ln-r*(un - g));
145 friction_law(e, rho);
146 auxN = lt - zt; ball_projection(auxN, rho);
147 for (i=0; i<N; ++i) t[i] = auxN[i] - e*no[i];
150 e = gmm::neg(-r*(un - g));
151 friction_law(e, rho);
152 auxN = - zt; ball_projection(auxN, rho);
153 for (i=0; i<N; ++i) t[i] = auxN[i] - e*no[i];
156 auxN = lnt - (r*(un-g) - f_coeff * gmm::vect_norm2(zt)) * no - zt;
157 De_Saxce_projection(auxN, no, f_coeff);
158 for (i=0; i<N; ++i) t[i] = auxN[i];
161 for (i=0; i<N; ++i) t[i] = lnt[i];
167 friction_law(e, rho);
168 auxN = lt; ball_projection(auxN, rho);
171 for (i=0; i<N; ++i) t[i] = auxN[i] - e*no[i];
174 auxN = lnt; De_Saxce_projection(auxN, no, f_coeff);
175 for (i=0; i<N; ++i) t[i] = auxN[i];
178 e = gmm::neg(ln-r*(un-g));
179 friction_law(e, rho);
180 auxN = zt - lt; ball_projection(auxN, rho); auxN += lt;
181 for (i=0; i<N; ++i) t[i] = ((e+ln)*no[i] + auxN[i])/ r;
184 e = r*(un-g) + gmm::pos(ln);
185 friction_law(gmm::neg(ln), rho);
186 auxN = lt; ball_projection(auxN, rho);
187 for (i=0; i<N; ++i) t[i] = (no[i]*e + zt[i] + lt[i] - auxN[i])/r;
191 De_Saxce_projection(auxN, no, f_coeff);
192 auxN -= lnt + (r*(un-g) - f_coeff * gmm::vect_norm2(zt)) * no + zt;
193 for (i=0; i<N; ++i) t[i] = -auxN[i]/r;
196 for (i=0; i<N; ++i) t[i] = -no[i];
200 for (i=0; i<N; ++i) t[i] = e*no[i];
203 e = -Heav(r*(un-g)-ln);
204 for (i=0; i<N; ++i) t[i] = e*no[i];
206 case UZAWA_PROJ_FRICT:
207 e = gmm::neg(ln - r*(un - g));
208 friction_law(e, rho);
209 auxN = lt - zt; ball_projection(auxN, rho);
210 for (i=0; i<N; ++i) t[i] = auxN[i] - e*no[i];
212 case UZAWA_PROJ_FRICT_SAXCE:
213 auxN = lnt - (r*(un-g) - f_coeff * gmm::vect_norm2(zt)) * no - zt;
214 De_Saxce_projection(auxN, no, f_coeff);
215 for (i=0; i<N; ++i) t[i] = auxN[i];
221 e = r * Heav(un - g);
222 for (i=0; i<N; ++i)
for (j=0; j<N; ++j) t[i*N+j] = e * no[i] * no[j];
225 e = r * Heav(r*(un - g)-ln);
226 for (i=0; i<N; ++i)
for (j=0; j<N; ++j) t[i*N+j] = e * no[i] * no[j];
230 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
231 t[i*N+j] = ((i == j) ? -scalar_type(1) : scalar_type(0));
234 friction_law(gmm::neg(ln), rho, rho_grad);
235 ball_projection_grad(lt, rho, GP);
237 coulomb = (rho_grad > 0) &&
bool(Heav(-ln));
238 if (coulomb) ball_projection_grad_r(lt, rho, V);
239 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
240 t[i*N+j] = no[i]*no[j]*e - GP(i,j) +
241 (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0));
244 augm_ln = ln - r*(un-g);
245 friction_law(gmm::neg(augm_ln), rho, rho_grad);
247 ball_projection_grad(auxN, rho, GP);
249 coulomb = (rho_grad > 0) &&
bool(Heav(-augm_ln));
250 if (coulomb) ball_projection_grad_r(auxN, rho, V);
251 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
252 t[i*N+j] = no[i]*no[j]*e - GP(i,j) +
253 (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0));
256 augm_ln = ln - r*(un-g);
257 friction_law(gmm::neg(augm_ln), rho, rho_grad);
259 ball_projection_grad(auxN, rho, GP); gmm::scale(GP, alpha);
261 coulomb = (rho_grad > 0) &&
bool(Heav(-augm_ln));
262 if (coulomb) ball_projection_grad_r(auxN, rho, V);
263 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
264 t[i*N+j] = no[i]*no[j]*e - GP(i,j) +
265 (coulomb ? rho_grad*V[i]*no[j] : scalar_type(0));
268 e =
alpha - scalar_type(1);
269 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
270 t[i*N+j] = no[i]*no[j]*e - ((i == j) ? alpha : scalar_type(0));
273 De_Saxce_projection_grad(lnt, no, f_coeff, GP);
274 for (i=0; i<N; ++i)
for (j=0; j<N; ++j) t[i*N+j] = -GP(j,i);
279 gmm::copy(gmm::identity_matrix(), GP); gmm::scale(GP, alpha);
280 gmm::rank_one_update(GP, gmm::scaled(no, scalar_type(1)-alpha), no);
281 if (nzt != scalar_type(0))
282 gmm::rank_one_update(GP, gmm::scaled(no, -f_coeff*alpha/nzt), zt);
283 for (i=0; i<N; ++i)
for (j=0; j<N; ++j) t[i*N+j] = - GP(i,j);
287 augm_ln = ln - r*(un-g);
288 friction_law(gmm::neg(augm_ln), rho, rho_grad);
290 ball_projection_grad(auxN, rho, GP);
292 coulomb = (rho_grad > 0) &&
bool(Heav(-augm_ln));
293 if (coulomb) ball_projection_grad_r(auxN, rho, V);
294 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
295 t[i*N+j] = (no[i]*no[j]*e + GP(i,j)
296 - ((i == j) ? scalar_type(1) : scalar_type(0))
297 - (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0))) / r;
300 friction_law(gmm::neg(ln), rho, rho_grad);
301 ball_projection_grad(lt, rho, GP);
303 coulomb = (rho_grad > 0) &&
bool(Heav(-ln));
304 if (coulomb) ball_projection_grad_r(lt, rho, V);
305 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
306 t[i*N+j] = (no[i]*no[j]*e + GP(i,j)
307 - ((i == j) ? scalar_type(1) : scalar_type(0))
308 - (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0))) / r;
311 De_Saxce_projection_grad(lnt, no, f_coeff, GP);
312 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
313 t[i*N+j] = (GP(i,j) - ((i == j) ? scalar_type(1) : scalar_type(0)))/r;
316 e = r * Heav(r*(un-g)-ln);
317 for (i=0; i<N; ++i)
for (j=0; j<N; ++j) t[i*N+j] = no[i]*no[j]*e;
320 friction_law(-ln, rho, rho_grad);
322 ball_projection_grad(auxN, rho, GP); gmm::scale(GP, alpha);
324 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
325 t[i*N+j] = r*(no[i]*no[j]*e + GP(i,j));
329 augm_ln = (option == K_UU_FRICT_V3) ? ln - r*(un-g) : - r*(un-g);
330 auxN = (option == K_UU_FRICT_V3) ? lt - zt : -zt;
331 friction_law(gmm::neg(augm_ln), rho, rho_grad);
332 ball_projection_grad(auxN, rho, GP); gmm::scale(GP, alpha);
334 coulomb = (rho_grad > 0) &&
bool(Heav(-augm_ln));
335 if (coulomb) ball_projection_grad_r(auxN, rho, V);
336 for (i=0; i<N; ++i)
for (j=0; j<N; ++j)
337 t[i*N+j] = r*(no[i]*no[j]*e + GP(i,j)
338 - (coulomb ? rho_grad*no[i]*V[j] : scalar_type(0)));
343 auxN = lnt - (r*(un-g) - f_coeff * nzt) * no - zt;
344 base_matrix A(N, N), B(N, N);
345 De_Saxce_projection_grad(auxN, no, f_coeff, A);
346 gmm::copy(gmm::identity_matrix(), B); gmm::scale(B, alpha);
347 gmm::rank_one_update(B, gmm::scaled(no, scalar_type(1)-alpha), no);
348 if (nzt != scalar_type(0))
349 gmm::rank_one_update(B, gmm::scaled(no, -f_coeff*alpha/nzt), zt);
351 for (i=0; i<N; ++i)
for (j=0; j<N; ++j) t[i*N+j] = r*GP(j,i);
354 default : GMM_ASSERT1(
false,
"Invalid option");
365 void contact_rigid_obstacle_nonlinear_term::prepare
366 (fem_interpolation_context& ctx,
size_type nb) {
372 ctx.pf()->interpolation(ctx, coeff, V, N);
375 if (gmm::vect_size(WT) == gmm::vect_size(U)) {
377 ctx.pf()->interpolation(ctx, coeff, auxN, N);
379 if (gmm::vect_size(VT) == gmm::vect_size(U)) {
381 ctx.pf()->interpolation(ctx, coeff, vt, N);
384 zt = (((V - un * no) - auxN) *
alpha + vt * (1 - gamma)) * r;
387 zt = ((V - un * no) - auxN) * (r *
alpha);
391 zt = (V - un * no) * (r * alpha);
399 ctx.pf()->interpolation_grad(ctx, coeff, grad, 1);
402 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
405 if (!contact_only && pmf_lambda) {
416 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
419 ctx.pf()->interpolation(ctx, coeff, lnt, N);
427 GMM_ASSERT1(!contact_only,
"Invalid friction option");
430 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
432 if (gmm::vect_size(tau_adhesion)) {
434 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
436 if (gmm::vect_size(tresca_limit)) {
438 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
439 tresca_lim = aux1[0];
445 default : GMM_ASSERT1(
false,
"Invalid option");
457 void contact_nonmatching_meshes_nonlinear_term::prepare
458 (fem_interpolation_context& ctx,
size_type nb) {
474 ctx.pf()->interpolation(ctx, coeff, V, N);
478 if (gmm::vect_size(WT1) == gmm::vect_size(U1)) {
480 ctx.pf()->interpolation(ctx, coeff, auxN, N);
482 zt = ((V - un1 * no) - auxN) * (r *
alpha) - zt;
484 zt = (V - un1 * no) * (r * alpha) - zt;
495 const projected_fem &pfe =
dynamic_cast<const projected_fem&
>(*ctx.pf());
496 pfe.projection_data(ctx, no, g);
497 gmm::scale(no, scalar_type(-1));
500 if (!contact_only && pmf_lambda) {
506 ctx.pf()->interpolation(ctx, coeff, V, N);
509 if (gmm::vect_size(WT2) == gmm::vect_size(U2)) {
511 ctx.pf()->interpolation(ctx, coeff, auxN, N);
513 zt = ((V - un * no) - auxN) * (r *
alpha);
515 zt = (V - un * no) * (r * alpha);
524 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
527 ctx.pf()->interpolation(ctx, coeff, lnt, N);
535 GMM_ASSERT1(!contact_only,
"Invalid friction option");
538 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
540 if (gmm::vect_size(tau_adhesion)) {
542 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
544 if (gmm::vect_size(tresca_limit)) {
546 ctx.pf()->interpolation(ctx, coeff, aux1, 1);
547 tresca_lim = aux1[0];
553 default : GMM_ASSERT1(
false,
"Invalid option");
565 template<
typename MAT,
typename VECT1>
566 void asm_Alart_Curnier_contact_rigid_obstacle_tangent_matrix
567 (MAT &Kul, MAT &Klu, MAT &Kll, MAT &Kuu,
572 scalar_type r,
const mesh_region &rg,
int option = 1) {
574 size_type subterm1 = (option == 3) ? K_UL_V2 : K_UL_V1;
575 size_type subterm2 = (option == 3) ? K_UL_V1 : K_UL_V3;
576 size_type subterm3 = (option == 3) ? K_LL_V2 : K_LL_V1;
577 size_type subterm4 = (option == 2) ? K_UU_V2 : K_UU_V1;
579 contact_rigid_obstacle_nonlinear_term
580 nterm1(subterm1, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda),
581 nterm2(subterm2, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda),
582 nterm3(subterm3, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda),
583 nterm4(subterm4, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda);
589 (
"M$1(#1,#3)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1).Base(#3))(i,:,i,:); "
590 "M$2(#3,#1)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3).vBase(#1))(i,:,:,i); "
591 "M$3(#3,#3)+=comp(NonLin$3(#1,#1,#2,#3).Base(#3).Base(#3))(i,:,:)");
595 (
"M$1(#1,#3)+=comp(NonLin$2(#1,#1,#2,#3).vBase(#1).Base(#3))(i,:,i,:); "
596 "M$3(#3,#3)+=comp(NonLin$3(#1,#1,#2,#3).Base(#3).Base(#3))(i,:,:);"
597 "M$4(#1,#1)+=comp(NonLin$4(#1,#1,#2,#3).vBase(#1).vBase(#1))(i,j,:,i,:,j)");
615 template<
typename MAT,
typename VECT1>
616 void asm_Alart_Curnier_contact_rigid_obstacle_tangent_matrix
617 (MAT &Kul, MAT &Klu, MAT &Kll, MAT &Kuu,
623 scalar_type
alpha,
const VECT1 *WT,
624 scalar_type gamma,
const VECT1 *VT,
625 const mesh_region &rg,
int option = 1) {
629 case 1 : subterm1 = K_UL_FRICT_V1; subterm2 = K_UL_FRICT_V4;
630 subterm3 = K_LL_FRICT_V1;
break;
631 case 2 : subterm1 = K_UL_FRICT_V3; subterm2 = K_UL_FRICT_V4;
632 subterm3 = K_LL_FRICT_V1;
break;
633 case 3 : subterm1 = K_UL_FRICT_V2; subterm2 = K_UL_FRICT_V5;
634 subterm3 = K_LL_FRICT_V2;
break;
635 case 4 : subterm1 = K_UL_FRICT_V7; subterm2 = K_UL_FRICT_V8;
636 subterm3 = K_LL_FRICT_V4;
break;
637 default : GMM_ASSERT1(
false,
"Incorrect option");
642 contact_rigid_obstacle_nonlinear_term
643 nterm1(subterm1, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
644 pmf_coeff, f_coeffs, alpha, WT, gamma, VT),
645 nterm2(subterm2, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
646 pmf_coeff, f_coeffs, alpha, WT, gamma, VT),
647 nterm3(subterm3, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
648 pmf_coeff, f_coeffs, alpha, WT, gamma, VT),
649 nterm4(subterm4, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
650 pmf_coeff, f_coeffs, alpha, WT, gamma, VT);
652 const std::string aux_fems = pmf_coeff ?
"#1,#2,#3,#4" :
"#1,#2,#3";
656 case 1:
case 3:
case 4:
658 (
"M$1(#1,#3)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1).vBase(#3))(i,j,:,i,:,j); "
659 "M$2(#3,#1)+=comp(NonLin$2(#1," + aux_fems +
").vBase(#3).vBase(#1))(i,j,:,j,:,i); "
660 "M$3(#3,#3)+=comp(NonLin$3(#1," + aux_fems +
").vBase(#3).vBase(#3))(i,j,:,i,:,j)");
664 (
"M$1(#1,#3)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1).vBase(#3))(i,j,:,i,:,j); "
665 "M$2(#3,#1)+=comp(NonLin$2(#1," + aux_fems +
").vBase(#3).vBase(#1))(i,j,:,j,:,i); "
666 "M$3(#3,#3)+=comp(NonLin$3(#1," + aux_fems +
").vBase(#3).vBase(#3))(i,j,:,i,:,j);"
667 "M$4(#1,#1)+=comp(NonLin$4(#1," + aux_fems +
").vBase(#1).vBase(#1))(i,j,:,i,:,j)");
687 template<
typename VECT1>
688 void asm_Alart_Curnier_contact_rigid_obstacle_rhs
689 (VECT1 &Ru, VECT1 &Rl,
694 scalar_type r,
const mesh_region &rg,
int option = 1) {
698 case 1 : subterm1 = RHS_U_V1;
break;
699 case 2 : subterm1 = RHS_U_V2;
break;
700 case 3 : subterm1 = RHS_U_V4;
break;
701 default : GMM_ASSERT1(
false,
"Incorrect option");
703 size_type subterm2 = (option == 3) ? RHS_L_V2 : RHS_L_V1;
705 contact_rigid_obstacle_nonlinear_term
706 nterm1(subterm1, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda),
707 nterm2(subterm2, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda);
710 assem.set(
"V$1(#1)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1))(i,:,i); "
711 "V$2(#3)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3))(i,:)");
724 template<
typename VECT1>
725 void asm_Alart_Curnier_contact_rigid_obstacle_rhs
726 (VECT1 &Ru, VECT1 &Rl,
732 scalar_type
alpha,
const VECT1 *WT,
733 scalar_type gamma,
const VECT1 *VT,
734 const mesh_region &rg,
int option = 1) {
738 case 1 : subterm1 = RHS_U_FRICT_V1; subterm2 = RHS_L_FRICT_V1;
break;
739 case 2 : subterm1 = RHS_U_FRICT_V6; subterm2 = RHS_L_FRICT_V1;
break;
740 case 3 : subterm1 = RHS_U_FRICT_V4; subterm2 = RHS_L_FRICT_V2;
break;
741 case 4 : subterm1 = RHS_U_FRICT_V5; subterm2 = RHS_L_FRICT_V4;
break;
742 default : GMM_ASSERT1(
false,
"Incorrect option");
745 contact_rigid_obstacle_nonlinear_term
746 nterm1(subterm1, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
747 pmf_coeff, f_coeffs, alpha, WT, gamma, VT),
748 nterm2(subterm2, r, mf_u, U, mf_obs, obs, &mf_lambda, &lambda,
749 pmf_coeff, f_coeffs, alpha, WT, gamma, VT);
751 const std::string aux_fems = pmf_coeff ?
"#1,#2,#3,#4" :
"#1,#2,#3";
754 assem.set(
"V$1(#1)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1))(i,:,i); "
755 "V$2(#3)+=comp(NonLin$2(#1," + aux_fems +
").vBase(#3))(i,:,i)");
769 struct integral_contact_rigid_obstacle_brick :
public virtual_brick {
780 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
781 const model::varnamelist &vl,
782 const model::varnamelist &dl,
783 const model::mimlist &mims,
784 model::real_matlist &matl,
785 model::real_veclist &vecl,
786 model::real_veclist &,
788 build_version version)
const {
789 GMM_ASSERT1(mims.size() == 1,
790 "Integral contact with rigid obstacle bricks need a single mesh_im");
791 GMM_ASSERT1(vl.size() == 2,
792 "Integral contact with rigid obstacle bricks need two variables");
793 GMM_ASSERT1(dl.size() >= 2 && dl.size() <= 7,
794 "Wrong number of data for integral contact with rigid obstacle "
795 <<
"brick, " << dl.size() <<
" should be between 2 and 7.");
796 GMM_ASSERT1(matl.size() ==
size_type(3 + (option == 2 && !contact_only)),
797 "Wrong number of terms for "
798 "integral contact with rigid obstacle brick");
809 const model_real_plain_vector &u = md.real_variable(vl[0]);
810 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
811 const model_real_plain_vector &lambda = md.real_variable(vl[1]);
812 const mesh_fem &mf_lambda = md.mesh_fem_of_variable(vl[1]);
813 GMM_ASSERT1(mf_lambda.get_qdim() == (contact_only ? 1 : mf_u.get_qdim()),
814 "The contact stress has not the right dimension");
815 const model_real_plain_vector &obstacle = md.real_variable(dl[0]);
816 const mesh_fem &mf_obstacle = md.mesh_fem_of_variable(dl[0]);
817 size_type sl = gmm::vect_size(obstacle) * mf_obstacle.get_qdim()
818 / mf_obstacle.nb_dof();
819 GMM_ASSERT1(sl == 1,
"the data corresponding to the obstacle has not "
822 const model_real_plain_vector &vr = md.real_variable(dl[1]);
823 GMM_ASSERT1(gmm::vect_size(vr) == 1,
"Parameter r should be a scalar");
824 const mesh_im &mim = *mims[0];
826 const model_real_plain_vector dummy_vec(0);
827 const model_real_plain_vector &friction_coeffs = contact_only
828 ? dummy_vec : md.real_variable(dl[2]);
829 const mesh_fem *pmf_coeff = contact_only ? 0 : md.pmesh_fem_of_variable(dl[2]);
830 sl = gmm::vect_size(friction_coeffs);
831 if (pmf_coeff) { sl *= pmf_coeff->get_qdim(); sl /= pmf_coeff->nb_dof(); }
832 GMM_ASSERT1(sl == 1 || sl == 2 || sl == 3 || contact_only,
833 "the data corresponding to the friction coefficient "
834 "has not the right format");
836 scalar_type
alpha = 1;
837 if (!contact_only && dl.size() >= 4) {
838 alpha = md.real_variable(dl[3])[0];
839 GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[3])) == 1,
840 "Parameter alpha should be a scalar");
843 const model_real_plain_vector *WT = 0;
844 if (!contact_only && dl.size() >= 5) {
845 if (dl[4].compare(vl[0]) != 0)
846 WT = &(md.real_variable(dl[4]));
847 else if (md.n_iter_of_variable(vl[0]) > 1)
848 WT = &(md.real_variable(vl[0],1));
851 scalar_type gamma = 1;
852 if (!contact_only && dl.size() >= 6) {
853 GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[5])) == 1,
854 "Parameter gamma should be a scalar");
855 gamma = md.real_variable(dl[5])[0];
858 const model_real_plain_vector *VT
859 = (!contact_only && dl.size()>=7) ? &(md.real_variable(dl[6])) : 0;
861 mesh_region rg(region);
862 mf_u.linked_mesh().intersect_with_mpi_region(rg);
864 if (version & model::BUILD_MATRIX) {
865 GMM_TRACE2(
"Integral contact with rigid obstacle friction tangent term");
868 size_type fourthmat = (matl.size() >= 4) ? 3 : 1;
870 asm_Alart_Curnier_contact_rigid_obstacle_tangent_matrix
871 (matl[0], matl[1], matl[2], matl[fourthmat], mim,
872 mf_u, u, mf_obstacle, obstacle, mf_lambda, lambda, vr[0],
875 asm_Alart_Curnier_contact_rigid_obstacle_tangent_matrix
876 (matl[0], matl[1], matl[2], matl[fourthmat], mim,
877 mf_u, u, mf_obstacle, obstacle, mf_lambda, lambda,
878 pmf_coeff, &friction_coeffs, vr[0], alpha, WT, gamma, VT,
882 if (version & model::BUILD_RHS) {
887 asm_Alart_Curnier_contact_rigid_obstacle_rhs
888 (vecl[0], vecl[2], mim,
889 mf_u, u, mf_obstacle, obstacle, mf_lambda, lambda, vr[0],
892 asm_Alart_Curnier_contact_rigid_obstacle_rhs
893 (vecl[0], vecl[2], mim,
894 mf_u, u, mf_obstacle, obstacle, mf_lambda, lambda,
895 pmf_coeff, &friction_coeffs, vr[0], alpha, WT, gamma, VT,
901 integral_contact_rigid_obstacle_brick(
bool contact_only_,
int option_) {
903 contact_only = contact_only_;
904 set_flags(contact_only
905 ?
"Integral contact with rigid obstacle brick"
906 :
"Integral contact and friction with rigid obstacle brick",
908 (option==2) && contact_only ,
922 (
model &md,
const mesh_im &mim,
const std::string &varname_u,
923 const std::string &multname_n,
const std::string &dataname_obs,
924 const std::string &dataname_r,
size_type region,
int option) {
927 = std::make_shared<integral_contact_rigid_obstacle_brick>(
true, option);
933 tl.push_back(model::term_description(varname_u, multname_n,
false));
934 tl.push_back(model::term_description(multname_n, varname_u,
false));
935 tl.push_back(model::term_description(multname_n, multname_n,
true));
938 tl.push_back(model::term_description(varname_u, multname_n,
true));
939 tl.push_back(model::term_description(varname_u, varname_u,
true));
940 tl.push_back(model::term_description(multname_n, multname_n,
true));
942 default :GMM_ASSERT1(
false,
943 "Incorrect option for integral contact brick");
946 model::varnamelist dl(1, dataname_obs);
947 dl.push_back(dataname_r);
949 model::varnamelist vl(1, varname_u);
950 vl.push_back(multname_n);
952 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
962 (
model &md,
const mesh_im &mim,
const std::string &varname_u,
963 const std::string &multname,
const std::string &dataname_obs,
964 const std::string &dataname_r,
const std::string &dataname_friction_coeffs,
966 const std::string &dataname_alpha,
const std::string &dataname_wt,
967 const std::string &dataname_gamma,
const std::string &dataname_vt) {
969 pbrick pbr = std::make_shared<integral_contact_rigid_obstacle_brick>
975 case 1:
case 3:
case 4:
976 tl.push_back(model::term_description(varname_u, multname,
false));
977 tl.push_back(model::term_description(multname, varname_u,
false));
978 tl.push_back(model::term_description(multname, multname,
true));
981 tl.push_back(model::term_description(varname_u, multname,
false));
982 tl.push_back(model::term_description(multname, varname_u,
false));
983 tl.push_back(model::term_description(multname, multname,
true));
984 tl.push_back(model::term_description(varname_u, varname_u,
true));
986 default :GMM_ASSERT1(
false,
987 "Incorrect option for integral contact brick");
989 model::varnamelist dl(1, dataname_obs);
990 dl.push_back(dataname_r);
991 dl.push_back(dataname_friction_coeffs);
992 if (dataname_alpha.size()) {
993 dl.push_back(dataname_alpha);
994 if (dataname_wt.size()) {
995 dl.push_back(dataname_wt);
996 if (dataname_gamma.size()) {
997 dl.push_back(dataname_gamma);
998 if (dataname_vt.size()) dl.push_back(dataname_vt);
1003 model::varnamelist vl(1, varname_u);
1004 vl.push_back(multname);
1006 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1016 template<
typename MAT,
typename VECT1>
1017 void asm_penalized_contact_rigid_obstacle_tangent_matrix
1023 scalar_type r,
const mesh_region &rg,
int option = 1) {
1025 contact_rigid_obstacle_nonlinear_term
1026 nterm((option == 1) ? K_UU_V1 : K_UU_V2, r,
1027 mf_u, U, mf_obs, obs, pmf_lambda, lambda);
1029 const std::string aux_fems = pmf_lambda ?
"#1,#2,#3":
"#1,#2";
1031 assem.set(
"M(#1,#1)+=comp(NonLin(#1," + aux_fems +
").vBase(#1).vBase(#1))(i,j,:,i,:,j)");
1043 template<
typename VECT1>
1044 void asm_penalized_contact_rigid_obstacle_rhs
1050 scalar_type r,
const mesh_region &rg,
int option = 1) {
1052 contact_rigid_obstacle_nonlinear_term
1053 nterm((option == 1) ? RHS_U_V5 : RHS_U_V2, r,
1054 mf_u, U, mf_obs, obs, pmf_lambda, lambda);
1056 const std::string aux_fems = pmf_lambda ?
"#1,#2,#3":
"#1,#2";
1058 assem.set(
"V(#1)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1))(i,:,i); ");
1069 template<
typename MAT,
typename VECT1>
1070 void asm_penalized_contact_rigid_obstacle_tangent_matrix
1077 scalar_type
alpha,
const VECT1 *WT,
1078 const mesh_region &rg,
int option = 1) {
1082 case 1 : subterm = K_UU_FRICT_V4;
break;
1083 case 2 : subterm = K_UU_FRICT_V3;
break;
1084 case 3 : subterm = K_UU_FRICT_V5;
break;
1087 contact_rigid_obstacle_nonlinear_term
1088 nterm(subterm, r, mf_u, U, mf_obs, obs, pmf_lambda, lambda,
1089 pmf_coeff, f_coeffs, alpha, WT);
1091 const std::string aux_fems = pmf_coeff ?
"#1,#2,#3,#4"
1092 : (pmf_lambda ?
"#1,#2,#3":
"#1,#2");
1094 assem.set(
"M(#1,#1)+=comp(NonLin(#1," + aux_fems +
").vBase(#1).vBase(#1))(i,j,:,i,:,j)");
1113 template<
typename VECT1>
1114 void asm_penalized_contact_rigid_obstacle_rhs
1121 scalar_type
alpha,
const VECT1 *WT,
1122 const mesh_region &rg,
int option = 1) {
1126 case 1 : subterm = RHS_U_FRICT_V7;
break;
1127 case 2 : subterm = RHS_U_FRICT_V6;
break;
1128 case 3 : subterm = RHS_U_FRICT_V8;
break;
1131 contact_rigid_obstacle_nonlinear_term
1132 nterm(subterm, r, mf_u, U, mf_obs, obs, pmf_lambda, lambda,
1133 pmf_coeff, f_coeffs, alpha, WT);
1135 const std::string aux_fems = pmf_coeff ?
"#1,#2,#3,#4"
1136 : (pmf_lambda ?
"#1,#2,#3":
"#1,#2");
1138 assem.set(
"V(#1)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1))(i,:,i); ");
1157 struct penalized_contact_rigid_obstacle_brick :
public virtual_brick {
1162 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
1163 const model::varnamelist &vl,
1164 const model::varnamelist &dl,
1165 const model::mimlist &mims,
1166 model::real_matlist &matl,
1167 model::real_veclist &vecl,
1168 model::real_veclist &,
1170 build_version version)
const {
1172 GMM_ASSERT1(mims.size() == 1,
1173 "Penalized contact with rigid obstacle bricks need a single mesh_im");
1174 const mesh_im &mim = *mims[0];
1177 GMM_ASSERT1(vl.size() == 1,
1178 "Penalized contact with rigid obstacle bricks need a single variable");
1179 const model_real_plain_vector &u = md.real_variable(vl[0]);
1180 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
1185 size_type nb_data_1 = ((option == 1) ? 2 : 3) + (contact_only ? 0 : 1);
1186 size_type nb_data_2 = nb_data_1 + (contact_only ? 0 : 2);
1187 GMM_ASSERT1(dl.size() >= nb_data_1 && dl.size() <= nb_data_2,
1188 "Wrong number of data for penalized contact with rigid obstacle "
1189 <<
"brick, " << dl.size() <<
" should be between "
1190 << nb_data_1 <<
" and " << nb_data_2 <<
".");
1193 const model_real_plain_vector &obs = md.real_variable(dl[nd]);
1194 const mesh_fem &mf_obs = md.mesh_fem_of_variable(dl[nd]);
1195 size_type sl = gmm::vect_size(obs) * mf_obs.get_qdim() / mf_obs.nb_dof();
1196 GMM_ASSERT1(sl == 1,
"the data corresponding to the obstacle has not "
1197 "the right format");
1200 const model_real_plain_vector &vr = md.real_variable(dl[nd]);
1201 GMM_ASSERT1(gmm::vect_size(vr) == 1,
"Parameter r should be a scalar");
1203 const model_real_plain_vector *lambda = 0;
1204 const mesh_fem *pmf_lambda = 0;
1207 lambda = &(md.real_variable(dl[nd]));
1208 pmf_lambda = md.pmesh_fem_of_variable(dl[nd]);
1209 sl = gmm::vect_size(*lambda) * pmf_lambda->get_qdim() / pmf_lambda->nb_dof();
1210 GMM_ASSERT1(sl == (contact_only ? 1 : N),
1211 "the data corresponding to the contact stress "
1212 "has not the right format");
1215 const model_real_plain_vector *f_coeffs = 0;
1216 const mesh_fem *pmf_coeff = 0;
1217 scalar_type
alpha = 1;
1218 const model_real_plain_vector *WT = 0;
1219 if (!contact_only) {
1221 f_coeffs = &(md.real_variable(dl[nd]));
1222 pmf_coeff = md.pmesh_fem_of_variable(dl[nd]);
1223 sl = gmm::vect_size(*f_coeffs);
1224 if (pmf_coeff) { sl *= pmf_coeff->get_qdim(); sl /= pmf_coeff->nb_dof(); }
1225 GMM_ASSERT1(sl == 1 || sl == 2 || sl == 3,
1226 "the data corresponding to the friction coefficient "
1227 "has not the right format");
1228 if (dl.size() > nd+1) {
1230 alpha = md.real_variable(dl[nd])[0];
1231 GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[nd])) == 1,
1232 "Parameter alpha should be a scalar");
1235 if (dl.size() > nd+1) {
1237 if (dl[nd].compare(vl[0]) != 0)
1238 WT = &(md.real_variable(dl[nd]));
1239 else if (md.n_iter_of_variable(vl[0]) > 1)
1240 WT = &(md.real_variable(vl[0],1));
1244 GMM_ASSERT1(matl.size() == 1,
"Wrong number of terms for "
1245 "penalized contact with rigid obstacle brick");
1247 mesh_region rg(region);
1248 mf_u.linked_mesh().intersect_with_mpi_region(rg);
1250 if (version & model::BUILD_MATRIX) {
1251 GMM_TRACE2(
"Penalized contact with rigid obstacle tangent term");
1254 asm_penalized_contact_rigid_obstacle_tangent_matrix
1255 (matl[0], mim, mf_u, u, mf_obs, obs, pmf_lambda, lambda,
1258 asm_penalized_contact_rigid_obstacle_tangent_matrix
1259 (matl[0], mim, mf_u, u, mf_obs, obs, pmf_lambda, lambda,
1260 pmf_coeff, f_coeffs, vr[0], alpha, WT, rg, option);
1263 if (version & model::BUILD_RHS) {
1266 asm_penalized_contact_rigid_obstacle_rhs
1267 (vecl[0], mim, mf_u, u, mf_obs, obs, pmf_lambda, lambda,
1270 asm_penalized_contact_rigid_obstacle_rhs
1271 (vecl[0], mim, mf_u, u, mf_obs, obs, pmf_lambda, lambda,
1272 pmf_coeff, f_coeffs, vr[0], alpha, WT, rg, option);
1277 penalized_contact_rigid_obstacle_brick(
bool contact_only_,
int option_) {
1278 contact_only = contact_only_;
1280 set_flags(contact_only
1281 ?
"Integral penalized contact with rigid obstacle brick"
1282 :
"Integral penalized contact and friction with rigid obstacle brick",
1283 false , contact_only ,
1297 (
model &md,
const mesh_im &mim,
const std::string &varname_u,
1298 const std::string &dataname_obs,
const std::string &dataname_r,
1299 size_type region,
int option,
const std::string &dataname_n) {
1301 pbrick pbr = std::make_shared<penalized_contact_rigid_obstacle_brick>
1305 tl.push_back(model::term_description(varname_u, varname_u,
true));
1307 model::varnamelist dl(1, dataname_obs);
1308 dl.push_back(dataname_r);
1311 case 2: dl.push_back(dataname_n);
break;
1312 default: GMM_ASSERT1(
false,
"Penalized contact brick : invalid option");
1315 model::varnamelist vl(1, varname_u);
1317 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1326 (
model &md,
const mesh_im &mim,
const std::string &varname_u,
1327 const std::string &dataname_obs,
const std::string &dataname_r,
1328 const std::string &dataname_friction_coeffs,
1329 size_type region,
int option,
const std::string &dataname_lambda,
1330 const std::string &dataname_alpha,
const std::string &dataname_wt) {
1332 pbrick pbr = std::make_shared<penalized_contact_rigid_obstacle_brick>
1336 tl.push_back(model::term_description(varname_u, varname_u,
false));
1338 model::varnamelist dl(1, dataname_obs);
1339 dl.push_back(dataname_r);
1342 case 2:
case 3: dl.push_back(dataname_lambda);
break;
1343 default: GMM_ASSERT1(
false,
"Penalized contact brick : invalid option");
1345 dl.push_back(dataname_friction_coeffs);
1346 if (dataname_alpha.size() > 0) {
1347 dl.push_back(dataname_alpha);
1348 if (dataname_wt.size() > 0) dl.push_back(dataname_wt);
1351 model::varnamelist vl(1, varname_u);
1353 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region);
1363 template<
typename MAT,
typename VEC>
1364 void asm_Alart_Curnier_contact_nonmatching_meshes_tangent_matrix
1365 (MAT &Ku1l, MAT &Klu1, MAT &Ku2l, MAT &Klu2, MAT &Kll,
1366 MAT &Ku1u1, MAT &Ku2u2, MAT &Ku1u2,
1371 scalar_type r,
const mesh_region &rg,
int option = 1) {
1373 size_type subterm1 = (option == 3) ? K_UL_V2 : K_UL_V1;
1374 size_type subterm2 = (option == 3) ? K_UL_V1 : K_UL_V3;
1375 size_type subterm3 = (option == 3) ? K_LL_V2 : K_LL_V1;
1376 size_type subterm4 = (option == 2) ? K_UU_V2 : K_UU_V1;
1378 contact_nonmatching_meshes_nonlinear_term
1379 nterm1(subterm1, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda),
1380 nterm2(subterm2, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda),
1381 nterm3(subterm3, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda),
1382 nterm4(subterm4, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda);
1388 (
"M$1(#1,#3)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1).Base(#3))(i,:,i,:); "
1389 "M$2(#3,#1)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3).vBase(#1))(i,:,:,i); "
1390 "M$3(#2,#3)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#2).Base(#3))(i,:,i,:); "
1391 "M$4(#3,#2)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3).vBase(#2))(i,:,:,i); "
1392 "M$5(#3,#3)+=comp(NonLin$3(#1,#1,#2,#3).Base(#3).Base(#3))(i,:,:)");
1396 (
"M$1(#1,#3)+=comp(NonLin$2(#1,#1,#2,#3).vBase(#1).Base(#3))(i,:,i,:); "
1397 "M$3(#2,#3)+=comp(NonLin$2(#1,#1,#2,#3).vBase(#2).Base(#3))(i,:,i,:); "
1398 "M$5(#3,#3)+=comp(NonLin$3(#1,#1,#2,#3).Base(#3).Base(#3))(i,:,:); "
1399 "M$6(#1,#1)+=comp(NonLin$4(#1,#1,#2,#3).vBase(#1).vBase(#1))(i,j,:,i,:,j); "
1400 "M$7(#2,#2)+=comp(NonLin$4(#1,#1,#2,#3).vBase(#2).vBase(#2))(i,j,:,i,:,j); "
1401 "M$8(#1,#2)+=comp(NonLin$4(#1,#1,#2,#3).vBase(#1).vBase(#2))(i,j,:,i,:,j)");
1422 gmm::scale(Ku2l, scalar_type(-1));
1424 gmm::scale(Klu2, scalar_type(-1));
1425 gmm::scale(Ku1u2, scalar_type(-1));
1428 template<
typename MAT,
typename VEC>
1429 void asm_Alart_Curnier_contact_nonmatching_meshes_tangent_matrix
1430 (MAT &Ku1l, MAT &Klu1, MAT &Ku2l, MAT &Klu2, MAT &Kll,
1431 MAT &Ku1u1, MAT &Ku2u2, MAT &Ku1u2, MAT &Ku2u1,
1437 scalar_type r, scalar_type
alpha,
1438 const VEC *WT1,
const VEC *WT2,
1439 const mesh_region &rg,
int option = 1) {
1444 subterm1 = K_UL_FRICT_V1; subterm2 = K_UL_FRICT_V4; subterm3 = K_LL_FRICT_V1;
1447 subterm1 = K_UL_FRICT_V3; subterm2 = K_UL_FRICT_V4; subterm3 = K_LL_FRICT_V1;
1450 subterm1 = K_UL_FRICT_V2; subterm2 = K_UL_FRICT_V5; subterm3 = K_LL_FRICT_V2;
1453 subterm1 = K_UL_FRICT_V7; subterm2 = K_UL_FRICT_V8; subterm3 = K_LL_FRICT_V4;
1455 default : GMM_ASSERT1(
false,
"Incorrect option");
1460 contact_nonmatching_meshes_nonlinear_term
1461 nterm1(subterm1, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1462 pmf_coeff, f_coeffs, alpha, WT1, WT2),
1463 nterm2(subterm2, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1464 pmf_coeff, f_coeffs, alpha, WT1, WT2),
1465 nterm3(subterm3, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1466 pmf_coeff, f_coeffs, alpha, WT1, WT2),
1467 nterm4(subterm4, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1468 pmf_coeff, f_coeffs, alpha, WT1, WT2);
1470 const std::string aux_fems = pmf_coeff ?
"#1,#2,#3,#4" :
"#1,#2,#3";
1474 case 1:
case 3:
case 4:
1476 (
"M$1(#1,#3)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1).vBase(#3))(i,j,:,i,:,j); "
1477 "M$2(#3,#1)+=comp(NonLin$2(#1," + aux_fems +
").vBase(#3).vBase(#1))(i,j,:,j,:,i); "
1478 "M$3(#2,#3)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#2).vBase(#3))(i,j,:,i,:,j); "
1479 "M$4(#3,#2)+=comp(NonLin$2(#1," + aux_fems +
").vBase(#3).vBase(#2))(i,j,:,j,:,i); "
1480 "M$5(#3,#3)+=comp(NonLin$3(#1," + aux_fems +
").vBase(#3).vBase(#3))(i,j,:,i,:,j)");
1484 (
"M$1(#1,#3)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1).vBase(#3))(i,j,:,i,:,j); "
1485 "M$2(#3,#1)+=comp(NonLin$2(#1," + aux_fems +
").vBase(#3).vBase(#1))(i,j,:,j,:,i); "
1486 "M$3(#2,#3)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#2).vBase(#3))(i,j,:,i,:,j); "
1487 "M$4(#3,#2)+=comp(NonLin$2(#1," + aux_fems +
").vBase(#3).vBase(#2))(i,j,:,j,:,i); "
1488 "M$5(#3,#3)+=comp(NonLin$3(#1," + aux_fems +
").vBase(#3).vBase(#3))(i,j,:,i,:,j); "
1489 "M$6(#1,#1)+=comp(NonLin$4(#1," + aux_fems +
").vBase(#1).vBase(#1))(i,j,:,i,:,j); "
1490 "M$7(#2,#2)+=comp(NonLin$4(#1," + aux_fems +
").vBase(#2).vBase(#2))(i,j,:,i,:,j); "
1491 "M$8(#1,#2)+=comp(NonLin$4(#1," + aux_fems +
").vBase(#1).vBase(#2))(i,j,:,i,:,j); "
1492 "M$9(#2,#1)+=comp(NonLin$4(#1," + aux_fems +
").vBase(#2).vBase(#1))(i,j,:,i,:,j)");
1516 gmm::scale(Ku2l, scalar_type(-1));
1517 gmm::scale(Klu2, scalar_type(-1));
1518 gmm::scale(Ku1u2, scalar_type(-1));
1521 template<
typename VECT1>
1522 void asm_Alart_Curnier_contact_nonmatching_meshes_rhs
1523 (VECT1 &Ru1, VECT1 &Ru2, VECT1 &Rl,
1528 scalar_type r,
const mesh_region &rg,
int option = 1) {
1532 case 1 : subterm1 = RHS_U_V1;
break;
1533 case 2 : subterm1 = RHS_U_V2;
break;
1534 case 3 : subterm1 = RHS_U_V4;
break;
1535 default : GMM_ASSERT1(
false,
"Incorrect option");
1537 size_type subterm2 = (option == 3) ? RHS_L_V2 : RHS_L_V1;
1539 contact_nonmatching_meshes_nonlinear_term
1540 nterm1(subterm1, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda),
1541 nterm2(subterm2, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda);
1544 assem.set(
"V$1(#1)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#1))(i,:,i); "
1545 "V$2(#2)+=comp(NonLin$1(#1,#1,#2,#3).vBase(#2))(i,:,i); "
1546 "V$3(#3)+=comp(NonLin$2(#1,#1,#2,#3).Base(#3))(i,:)");
1558 gmm::scale(Ru2, scalar_type(-1));
1561 template<
typename VECT1>
1562 void asm_Alart_Curnier_contact_nonmatching_meshes_rhs
1563 (VECT1 &Ru1, VECT1 &Ru2, VECT1 &Rl,
1569 scalar_type r, scalar_type
alpha,
1570 const VECT1 *WT1,
const VECT1 *WT2,
1571 const mesh_region &rg,
int option = 1) {
1575 case 1 : subterm1 = RHS_U_FRICT_V1; subterm2 = RHS_L_FRICT_V1;
break;
1576 case 2 : subterm1 = RHS_U_FRICT_V6; subterm2 = RHS_L_FRICT_V1;
break;
1577 case 3 : subterm1 = RHS_U_FRICT_V4; subterm2 = RHS_L_FRICT_V2;
break;
1578 case 4 : subterm1 = RHS_U_FRICT_V5; subterm2 = RHS_L_FRICT_V4;
break;
1579 default : GMM_ASSERT1(
false,
"Incorrect option");
1582 contact_nonmatching_meshes_nonlinear_term
1583 nterm1(subterm1, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1584 pmf_coeff, f_coeffs, alpha, WT1, WT2),
1585 nterm2(subterm2, r, mf_u1, U1, mf_u2, U2, &mf_lambda, &lambda,
1586 pmf_coeff, f_coeffs, alpha, WT1, WT2);
1588 const std::string aux_fems = pmf_coeff ?
"#1,#2,#3,#4" :
"#1,#2,#3";
1591 assem.set(
"V$1(#1)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1))(i,:,i); "
1592 "V$2(#2)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#2))(i,:,i); "
1593 "V$3(#3)+=comp(NonLin$2(#1," + aux_fems +
").vBase(#3))(i,:,i)");
1607 gmm::scale(Ru2, scalar_type(-1));
1610 struct integral_contact_nonmatching_meshes_brick :
public virtual_brick {
1623 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
1624 const model::varnamelist &vl,
1625 const model::varnamelist &dl,
1626 const model::mimlist &mims,
1627 model::real_matlist &matl,
1628 model::real_veclist &vecl,
1629 model::real_veclist &,
1631 build_version version)
const {
1633 GMM_ASSERT1(mims.size() == 1,
1634 "Integral contact between nonmatching meshes bricks need a single mesh_im");
1635 const mesh_im &mim = *mims[0];
1640 GMM_ASSERT1(vl.size() == 3,
1641 "Integral contact between nonmatching meshes bricks need three variables");
1642 const model_real_plain_vector &u1 = md.real_variable(vl[0]);
1643 const model_real_plain_vector &u2 = md.real_variable(vl[1]);
1644 const mesh_fem &mf_u1 = md.mesh_fem_of_variable(vl[0]);
1645 const mesh_fem &mf_u2 = md.mesh_fem_of_variable(vl[1]);
1646 const model_real_plain_vector &lambda = md.real_variable(vl[2]);
1647 const mesh_fem &mf_lambda = md.mesh_fem_of_variable(vl[2]);
1648 GMM_ASSERT1(mf_lambda.get_qdim() == (contact_only ? 1 : mf_u1.get_qdim()),
1649 "The contact stress variable has not the right dimension");
1654 GMM_ASSERT1(dl.size() == 1,
1655 "Wrong number of data for integral contact between nonmatching meshes "
1656 <<
"brick, the number of data should be equal to 1 .");
1659 GMM_ASSERT1(dl.size() >= 2 && dl.size() <= 5,
1660 "Wrong number of data for integral contact between nonmatching meshes "
1661 <<
"brick, it should be between 2 and 5 instead of " << dl.size() <<
" .");
1663 const model_real_plain_vector &vr = md.real_variable(dl[0]);
1664 GMM_ASSERT1(gmm::vect_size(vr) == 1,
"Parameter r should be a scalar");
1666 const model_real_plain_vector *f_coeffs = 0, *WT1 = 0, *WT2 = 0;
1667 const mesh_fem *pmf_coeff = 0;
1668 scalar_type
alpha = 1;
1669 if (!contact_only) {
1670 f_coeffs = &(md.real_variable(dl[1]));
1671 pmf_coeff = md.pmesh_fem_of_variable(dl[1]);
1673 size_type sl = gmm::vect_size(*f_coeffs);
1674 if (pmf_coeff) { sl *= pmf_coeff->get_qdim(); sl /= pmf_coeff->nb_dof(); }
1675 GMM_ASSERT1(sl == 1 || sl == 2 || sl ==3,
1676 "the data corresponding to the friction coefficient "
1677 "has not the right format");
1679 if (dl.size() >= 3) {
1680 alpha = md.real_variable(dl[2])[0];
1681 GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[2])) == 1,
1682 "Parameter alpha should be a scalar");
1685 if (dl.size() >= 4) {
1686 if (dl[3].compare(vl[0]) != 0)
1687 WT1 = &(md.real_variable(dl[3]));
1688 else if (md.n_iter_of_variable(vl[0]) > 1)
1689 WT1 = &(md.real_variable(vl[0],1));
1692 if (dl.size() >= 5) {
1693 if (dl[4].compare(vl[1]) != 0)
1694 WT2 = &(md.real_variable(dl[4]));
1695 else if (md.n_iter_of_variable(vl[1]) > 1)
1696 WT2 = &(md.real_variable(vl[1],1));
1701 GMM_ASSERT1(matl.size() ==
size_type(3 +
1702 2 * !is_symmetric() +
1704 1 * (option == 2 && !is_symmetric())),
1705 "Wrong number of terms for "
1706 "integral contact between nonmatching meshes brick");
1708 mesh_region rg(region);
1709 mf_u1.linked_mesh().intersect_with_mpi_region(rg);
1711 size_type N = mf_u1.linked_mesh().dim();
1718 mf_u2_proj.set_finite_element(mim.linked_mesh().convex_index(), pfem_proj);
1721 size_type nbdof_lambda = mf_lambda.nb_dof();
1723 size_type nbsub = mf_u2_proj.nb_basic_dof();
1725 std::vector<size_type> ind;
1726 mf_u2_proj.get_global_dof_index(ind);
1727 gmm::unsorted_sub_index SUBI(ind);
1729 gmm::csc_matrix<scalar_type> Esub(nbsub, nbdof2);
1730 if (mf_u2.is_reduced())
1731 gmm::copy(gmm::sub_matrix(mf_u2.extension_matrix(),
1732 SUBI, gmm::sub_interval(0, nbdof2)),
1735 model_real_plain_vector u2_proj(nbsub);
1736 if (mf_u2.is_reduced())
1739 gmm::copy(gmm::sub_vector(u2, SUBI), u2_proj);
1741 model_real_plain_vector WT2_proj(0);
1744 if (mf_u2.is_reduced())
1747 gmm::copy(gmm::sub_vector(*WT2, SUBI), WT2_proj);
1760 if (version & model::BUILD_MATRIX) {
1761 GMM_TRACE2(
"Integral contact between nonmatching meshes "
1763 for (
size_type i = 0; i < matl.size(); i++) gmm::clear(matl[i]);
1765 model_real_sparse_matrix dummy_mat(0, 0);
1766 model_real_sparse_matrix &Klu1 = (LU1 ==
size_type(-1)) ? dummy_mat : matl[LU1];
1767 model_real_sparse_matrix &Ku1u1 = (U1U1 ==
size_type(-1)) ? dummy_mat : matl[U1U1];
1769 model_real_sparse_matrix Ku2l(nbsub, nbdof_lambda);
1770 model_real_sparse_matrix Klu2(nbdof_lambda, nbsub);
1771 model_real_sparse_matrix Ku2u2(nbsub, nbsub);
1772 model_real_sparse_matrix Ku1u2(nbdof1, nbsub);
1773 model_real_sparse_matrix Ku2u1(nbsub, nbdof1);
1776 asm_Alart_Curnier_contact_nonmatching_meshes_tangent_matrix
1777 (matl[U1L], Klu1, Ku2l, Klu2, matl[LL], Ku1u1, Ku2u2, Ku1u2,
1778 mim, mf_u1, u1, mf_u2_proj, u2_proj, mf_lambda, lambda,
1781 asm_Alart_Curnier_contact_nonmatching_meshes_tangent_matrix
1782 (matl[U1L], Klu1, Ku2l, Klu2, matl[LL], Ku1u1, Ku2u2, Ku1u2, Ku2u1,
1783 mim, mf_u1, u1, mf_u2_proj, u2_proj, mf_lambda, lambda,
1784 pmf_coeff, f_coeffs, vr[0], alpha, WT1, &WT2_proj, rg, option);
1786 if (mf_u2.is_reduced()) {
1787 gmm::mult(gmm::transposed(Esub), Ku2l, matl[U2L]);
1790 model_real_sparse_matrix tmp(nbsub, nbdof2);
1792 gmm::mult(gmm::transposed(Esub), tmp, matl[U2U2]);
1798 gmm::copy(Ku2l, gmm::sub_matrix(matl[U2L], SUBI, gmm::sub_interval(0, nbdof_lambda)));
1800 gmm::copy(Klu2, gmm::sub_matrix(matl[LU2], gmm::sub_interval(0, nbdof_lambda), SUBI));
1802 gmm::copy(Ku2u2, gmm::sub_matrix(matl[U2U2], SUBI));
1804 gmm::copy(Ku1u2, gmm::sub_matrix(matl[U1U2], gmm::sub_interval(0, nbdof1), SUBI));
1806 gmm::copy(Ku2u1, gmm::sub_matrix(matl[U2U1], SUBI, gmm::sub_interval(0, nbdof1)));
1810 if (version & model::BUILD_RHS) {
1811 for (
size_type i = 0; i < matl.size(); i++) gmm::clear(vecl[i]);
1813 model_real_plain_vector Ru2(nbsub);
1816 asm_Alart_Curnier_contact_nonmatching_meshes_rhs
1817 (vecl[U1L], Ru2, vecl[LL],
1818 mim, mf_u1, u1, mf_u2_proj, u2_proj, mf_lambda, lambda,
1821 asm_Alart_Curnier_contact_nonmatching_meshes_rhs
1822 (vecl[U1L], Ru2, vecl[LL],
1823 mim, mf_u1, u1, mf_u2_proj, u2_proj, mf_lambda, lambda,
1824 pmf_coeff, f_coeffs, vr[0], alpha, WT1, &WT2_proj, rg, option);
1826 if (mf_u2.is_reduced())
1827 gmm::mult(gmm::transposed(Esub), Ru2, vecl[U2L]);
1829 gmm::copy(Ru2, gmm::sub_vector(vecl[U2L], SUBI));
1834 bool contact_only_,
int option_)
1835 : rg1(rg1_), rg2(rg2_), pfem_proj(0),
1836 contact_only(contact_only_), option(option_)
1838 set_flags(contact_only
1839 ?
"Integral contact between nonmatching meshes brick"
1840 :
"Integral contact and friction between nonmatching "
1843 (option==2) && contact_only ,
1848 ~integral_contact_nonmatching_meshes_brick()
1860 (
model &md,
const mesh_im &mim,
const std::string &varname_u1,
1861 const std::string &varname_u2,
const std::string &multname_n,
1862 const std::string &dataname_r,
1865 pbrick pbr = std::make_shared<integral_contact_nonmatching_meshes_brick>
1866 (region1, region2,
true , option);
1872 tl.push_back(model::term_description(varname_u1, multname_n,
false));
1873 tl.push_back(model::term_description(multname_n, varname_u1,
false));
1874 tl.push_back(model::term_description(varname_u2, multname_n,
false));
1875 tl.push_back(model::term_description(multname_n, varname_u2,
false));
1876 tl.push_back(model::term_description(multname_n, multname_n,
true));
1879 tl.push_back(model::term_description(varname_u1, multname_n,
true));
1880 tl.push_back(model::term_description(varname_u2, multname_n,
true));
1881 tl.push_back(model::term_description(multname_n, multname_n,
true));
1882 tl.push_back(model::term_description(varname_u1, varname_u1,
true));
1883 tl.push_back(model::term_description(varname_u2, varname_u2,
true));
1884 tl.push_back(model::term_description(varname_u1, varname_u2,
true));
1886 default : GMM_ASSERT1(
false,
1887 "Incorrect option for integral contact brick");
1890 model::varnamelist dl(1, dataname_r);
1892 model::varnamelist vl(1, varname_u1);
1893 vl.push_back(varname_u2);
1894 vl.push_back(multname_n);
1896 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region1);
1906 (
model &md,
const mesh_im &mim,
const std::string &varname_u1,
1907 const std::string &varname_u2,
const std::string &multname,
1908 const std::string &dataname_r,
const std::string &dataname_friction_coeffs,
1910 const std::string &dataname_alpha,
1911 const std::string &dataname_wt1,
const std::string &dataname_wt2) {
1913 pbrick pbr = std::make_shared<integral_contact_nonmatching_meshes_brick>
1914 (region1, region2,
false , option);
1919 case 1 :
case 3 :
case 4 :
1920 tl.push_back(model::term_description(varname_u1, multname,
false));
1921 tl.push_back(model::term_description(multname, varname_u1,
false));
1922 tl.push_back(model::term_description(varname_u2, multname,
false));
1923 tl.push_back(model::term_description(multname, varname_u2,
false));
1924 tl.push_back(model::term_description(multname, multname,
true));
1927 tl.push_back(model::term_description(varname_u1, multname,
false));
1928 tl.push_back(model::term_description(multname, varname_u1,
false));
1929 tl.push_back(model::term_description(varname_u2, multname,
false));
1930 tl.push_back(model::term_description(multname, varname_u2,
false));
1931 tl.push_back(model::term_description(multname, multname,
true));
1932 tl.push_back(model::term_description(varname_u1, varname_u1,
true));
1933 tl.push_back(model::term_description(varname_u2, varname_u2,
true));
1934 tl.push_back(model::term_description(varname_u1, varname_u2,
true));
1935 tl.push_back(model::term_description(varname_u2, varname_u1,
true));
1937 default : GMM_ASSERT1(
false,
1938 "Incorrect option for integral contact brick");
1941 model::varnamelist dl(1, dataname_r);
1942 dl.push_back(dataname_friction_coeffs);
1943 if (dataname_alpha.size()) {
1944 dl.push_back(dataname_alpha);
1945 if (dataname_wt1.size()) {
1946 dl.push_back(dataname_wt1);
1947 if (dataname_wt2.size()) {
1948 dl.push_back(dataname_wt2);
1954 model::varnamelist vl(1, varname_u1);
1955 vl.push_back(varname_u2);
1956 vl.push_back(multname);
1958 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region1);
1968 template<
typename MAT,
typename VECT1>
1969 void asm_penalized_contact_nonmatching_meshes_tangent_matrix
1970 (MAT &Ku1u1, MAT &Ku2u2, MAT &Ku1u2,
1975 const scalar_type r,
const mesh_region &rg,
int option = 1) {
1977 contact_nonmatching_meshes_nonlinear_term
1978 nterm((option == 1) ? K_UU_V1 : K_UU_V2, r,
1979 mf_u1, U1, mf_u2, U2, pmf_lambda, lambda);
1981 const std::string aux_fems = pmf_lambda ?
"#1,#2,#3" :
"#1,#2";
1984 assem.set(
"M$1(#1,#1)+=comp(NonLin(#1," + aux_fems +
").vBase(#1).vBase(#1))(i,j,:,i,:,j); "
1985 "M$2(#2,#2)+=comp(NonLin(#1," + aux_fems +
").vBase(#2).vBase(#2))(i,j,:,i,:,j); "
1986 "M$3(#1,#2)+=comp(NonLin(#1," + aux_fems +
").vBase(#1).vBase(#2))(i,j,:,i,:,j)");
1998 gmm::scale(Ku1u2, scalar_type(-1));
2001 template<
typename VECT1>
2002 void asm_penalized_contact_nonmatching_meshes_rhs
2003 (VECT1 &Ru1, VECT1 &Ru2,
2008 scalar_type r,
const mesh_region &rg,
int option = 1) {
2010 contact_nonmatching_meshes_nonlinear_term
2011 nterm((option == 1) ? RHS_U_V5 : RHS_U_V2, r,
2012 mf_u1, U1, mf_u2, U2, pmf_lambda, lambda);
2014 const std::string aux_fems = pmf_lambda ?
"#1,#2,#3":
"#1,#2";
2016 assem.set(
"V$1(#1)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1))(i,:,i); "
2017 "V$2(#2)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#2))(i,:,i)");
2028 gmm::scale(Ru2, scalar_type(-1));
2032 template<
typename MAT,
typename VECT1>
2033 void asm_penalized_contact_nonmatching_meshes_tangent_matrix
2034 (MAT &Ku1u1, MAT &Ku2u2, MAT &Ku1u2, MAT &Ku2u1,
2040 scalar_type
alpha,
const VECT1 *WT1,
const VECT1 *WT2,
2041 const mesh_region &rg,
int option = 1) {
2045 case 1 : subterm = K_UU_FRICT_V4;
break;
2046 case 2 : subterm = K_UU_FRICT_V3;
break;
2047 case 3 : subterm = K_UU_FRICT_V5;
break;
2050 contact_nonmatching_meshes_nonlinear_term
2051 nterm(subterm, r, mf_u1, U1, mf_u2, U2, pmf_lambda, lambda,
2052 pmf_coeff, f_coeffs, alpha, WT1, WT2);
2054 const std::string aux_fems = pmf_coeff ?
"#1,#2,#3,#4"
2055 : (pmf_lambda ?
"#1,#2,#3":
"#1,#2");
2058 assem.set(
"M$1(#1,#1)+=comp(NonLin(#1," + aux_fems +
").vBase(#1).vBase(#1))(i,j,:,i,:,j); "
2059 "M$2(#2,#2)+=comp(NonLin(#1," + aux_fems +
").vBase(#2).vBase(#2))(i,j,:,i,:,j); "
2060 "M$3(#1,#2)+=comp(NonLin(#1," + aux_fems +
").vBase(#1).vBase(#2))(i,j,:,i,:,j); "
2061 "M$4(#2,#1)+=comp(NonLin(#1," + aux_fems +
").vBase(#2).vBase(#1))(i,j,:,i,:,j)");
2081 gmm::scale(Ku1u2, scalar_type(-1));
2082 gmm::scale(Ku2u1, scalar_type(-1));
2086 template<
typename VECT1>
2087 void asm_penalized_contact_nonmatching_meshes_rhs
2088 (VECT1 &Ru1, VECT1 &Ru2,
2094 scalar_type
alpha,
const VECT1 *WT1,
const VECT1 *WT2,
2095 const mesh_region &rg,
int option = 1) {
2099 case 1 : subterm = RHS_U_FRICT_V7;
break;
2100 case 2 : subterm = RHS_U_FRICT_V6;
break;
2101 case 3 : subterm = RHS_U_FRICT_V8;
break;
2104 contact_nonmatching_meshes_nonlinear_term
2105 nterm(subterm, r, mf_u1, U1, mf_u2, U2, pmf_lambda, lambda,
2106 pmf_coeff, f_coeffs, alpha, WT1, WT2);
2108 const std::string aux_fems = pmf_coeff ?
"#1,#2,#3,#4"
2109 : (pmf_lambda ?
"#1,#2,#3":
"#1,#2");
2111 assem.set(
"V$1(#1)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#1))(i,:,i); "
2112 "V$2(#2)+=comp(NonLin$1(#1," + aux_fems +
").vBase(#2))(i,:,i)");
2131 gmm::scale(Ru2, scalar_type(-1));
2135 struct penalized_contact_nonmatching_meshes_brick :
public virtual_brick {
2143 virtual void asm_real_tangent_terms(
const model &md,
size_type ,
2144 const model::varnamelist &vl,
2145 const model::varnamelist &dl,
2146 const model::mimlist &mims,
2147 model::real_matlist &matl,
2148 model::real_veclist &vecl,
2149 model::real_veclist &,
2151 build_version version)
const {
2153 GMM_ASSERT1(mims.size() == 1,
2154 "Penalized contact between nonmatching meshes bricks need a single mesh_im");
2155 const mesh_im &mim = *mims[0];
2158 GMM_ASSERT1(vl.size() == 2,
2159 "Penalized contact between nonmatching meshes bricks need two variables");
2160 const model_real_plain_vector &u1 = md.real_variable(vl[0]);
2161 const model_real_plain_vector &u2 = md.real_variable(vl[1]);
2162 const mesh_fem &mf_u1 = md.mesh_fem_of_variable(vl[0]);
2163 const mesh_fem &mf_u2 = md.mesh_fem_of_variable(vl[1]);
2165 size_type N = mf_u1.linked_mesh().dim();
2168 size_type nb_data_1 = ((option == 1) ? 1 : 2) + (contact_only ? 0 : 1);
2169 size_type nb_data_2 = nb_data_1 + (contact_only ? 0 : 3);
2170 GMM_ASSERT1(dl.size() >= nb_data_1 && dl.size() <= nb_data_2,
2171 "Wrong number of data for penalized contact between nonmatching meshes "
2172 <<
"brick, " << dl.size() <<
" should be between "
2173 << nb_data_1 <<
" and " << nb_data_2 <<
".");
2176 const model_real_plain_vector &vr = md.real_variable(dl[nd]);
2177 GMM_ASSERT1(gmm::vect_size(vr) == 1,
"Parameter r should be a scalar");
2180 const model_real_plain_vector *lambda = 0;
2181 const mesh_fem *pmf_lambda = 0;
2184 lambda = &(md.real_variable(dl[nd]));
2185 pmf_lambda = md.pmesh_fem_of_variable(dl[nd]);
2186 sl = gmm::vect_size(*lambda) * pmf_lambda->get_qdim() / pmf_lambda->nb_dof();
2187 GMM_ASSERT1(sl == (contact_only ? 1 : N),
2188 "the data corresponding to the contact stress "
2189 "has not the right format");
2192 const model_real_plain_vector *f_coeffs = 0;
2193 const mesh_fem *pmf_coeff = 0;
2194 scalar_type
alpha = 1;
2195 const model_real_plain_vector *WT1 = 0;
2196 const model_real_plain_vector *WT2 = 0;
2197 if (!contact_only) {
2199 f_coeffs = &(md.real_variable(dl[nd]));
2200 pmf_coeff = md.pmesh_fem_of_variable(dl[nd]);
2201 sl = gmm::vect_size(*f_coeffs);
2202 if (pmf_coeff) { sl *= pmf_coeff->get_qdim(); sl /= pmf_coeff->nb_dof(); }
2203 GMM_ASSERT1(sl == 1 || sl == 2 || sl == 3,
2204 "the data corresponding to the friction coefficient "
2205 "has not the right format");
2207 if (dl.size() > nd+1) {
2209 alpha = md.real_variable(dl[nd])[0];
2210 GMM_ASSERT1(gmm::vect_size(md.real_variable(dl[nd])) == 1,
2211 "Parameter alpha should be a scalar");
2214 if (dl.size() > nd+1) {
2216 if (dl[nd].compare(vl[0]) != 0)
2217 WT1 = &(md.real_variable(dl[nd]));
2218 else if (md.n_iter_of_variable(vl[0]) > 1)
2219 WT1 = &(md.real_variable(vl[0],1));
2222 if (dl.size() > nd+1) {
2224 if (dl[nd].compare(vl[1]) != 0)
2225 WT2 = &(md.real_variable(dl[nd]));
2226 else if (md.n_iter_of_variable(vl[1]) > 1)
2227 WT2 = &(md.real_variable(vl[1],1));
2231 GMM_ASSERT1(matl.size() == (contact_only ? 3 : 4),
2232 "Wrong number of terms for penalized contact "
2233 "between nonmatching meshes brick");
2235 mesh_region rg(region);
2236 mf_u1.linked_mesh().intersect_with_mpi_region(rg);
2243 mf_u2_proj.set_finite_element(mim.linked_mesh().convex_index(), pfem_proj);
2249 std::vector<size_type> ind;
2250 mf_u2_proj.get_global_dof_index(ind);
2251 gmm::unsorted_sub_index SUBI(ind);
2253 gmm::csc_matrix<scalar_type> Esub(nbsub, nbdof2);
2254 if (mf_u2.is_reduced())
2255 gmm::copy(gmm::sub_matrix(mf_u2.extension_matrix(),
2256 SUBI, gmm::sub_interval(0, nbdof2)),
2259 model_real_plain_vector u2_proj(nbsub);
2260 if (mf_u2.is_reduced())
2263 gmm::copy(gmm::sub_vector(u2, SUBI), u2_proj);
2265 model_real_plain_vector WT2_proj(0);
2268 if (mf_u2.is_reduced())
2271 gmm::copy(gmm::sub_vector(*WT2, SUBI), WT2_proj);
2274 if (version & model::BUILD_MATRIX) {
2275 GMM_TRACE2(
"Penalized contact between nonmatching meshes tangent term");
2280 model_real_sparse_matrix Ku2u2(nbsub,nbsub);
2281 model_real_sparse_matrix Ku1u2(nbdof1,nbsub);
2284 asm_penalized_contact_nonmatching_meshes_tangent_matrix
2285 (matl[0], Ku2u2, Ku1u2, mim, mf_u1, u1, mf_u2_proj, u2_proj,
2286 pmf_lambda, lambda, vr[0], rg, option);
2290 model_real_sparse_matrix Ku2u1(nbsub,nbdof1);
2291 asm_penalized_contact_nonmatching_meshes_tangent_matrix
2292 (matl[0], Ku2u2, Ku1u2, Ku2u1, mim, mf_u1, u1, mf_u2_proj, u2_proj,
2293 pmf_lambda, lambda, pmf_coeff, f_coeffs, vr[0], alpha, WT1, &WT2_proj, rg, option);
2294 if (mf_u2.is_reduced())
2295 gmm::mult(gmm::transposed(Esub), Ku2u1, matl[3]);
2297 gmm::copy(Ku2u1, gmm::sub_matrix(matl[3], SUBI, gmm::sub_interval(0, nbdof1)));
2300 if (mf_u2.is_reduced()) {
2301 model_real_sparse_matrix tmp(nbsub, nbdof2);
2303 gmm::mult(gmm::transposed(Esub), tmp, matl[1]);
2307 gmm::copy(Ku2u2, gmm::sub_matrix(matl[1], SUBI));
2308 gmm::copy(Ku1u2, gmm::sub_matrix(matl[2], gmm::sub_interval(0, nbdof1), SUBI));
2312 if (version & model::BUILD_RHS) {
2316 model_real_plain_vector Ru2(nbsub);
2319 asm_penalized_contact_nonmatching_meshes_rhs
2320 (vecl[0], Ru2, mim, mf_u1, u1, mf_u2_proj, u2_proj, pmf_lambda, lambda,
2323 asm_penalized_contact_nonmatching_meshes_rhs
2324 (vecl[0], Ru2, mim, mf_u1, u1, mf_u2_proj, u2_proj, pmf_lambda, lambda,
2325 pmf_coeff, f_coeffs, vr[0], alpha, WT1, &WT2_proj, rg, option);
2327 if (mf_u2.is_reduced())
2328 gmm::mult(gmm::transposed(Esub), Ru2, vecl[1]);
2330 gmm::copy(Ru2, gmm::sub_vector(vecl[1], SUBI));
2335 bool contact_only_,
int option_)
2336 : rg1(rg1_), rg2(rg2_), pfem_proj(0),
2337 contact_only(contact_only_), option(option_)
2339 set_flags(contact_only
2340 ?
"Integral penalized contact between nonmatching meshes brick"
2341 :
"Integral penalized contact and friction between nonmatching "
2343 false , contact_only ,
2348 ~penalized_contact_nonmatching_meshes_brick()
2360 (
model &md,
const mesh_im &mim,
const std::string &varname_u1,
2361 const std::string &varname_u2,
const std::string &dataname_r,
2363 int option,
const std::string &dataname_n) {
2365 pbrick pbr = std::make_shared<penalized_contact_nonmatching_meshes_brick>
2366 (region1, region2,
true , option);
2368 tl.push_back(model::term_description(varname_u1, varname_u1,
true));
2369 tl.push_back(model::term_description(varname_u2, varname_u2,
true));
2370 tl.push_back(model::term_description(varname_u1, varname_u2,
true));
2372 model::varnamelist dl(1, dataname_r);
2375 case 2: dl.push_back(dataname_n);
break;
2376 default: GMM_ASSERT1(
false,
"Penalized contact brick : invalid option");
2379 model::varnamelist vl(1, varname_u1);
2380 vl.push_back(varname_u2);
2382 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region1);
2392 (
model &md,
const mesh_im &mim,
const std::string &varname_u1,
2393 const std::string &varname_u2,
const std::string &dataname_r,
2394 const std::string &dataname_friction_coeffs,
2396 const std::string &dataname_lambda,
const std::string &dataname_alpha,
2397 const std::string &dataname_wt1,
const std::string &dataname_wt2) {
2399 pbrick pbr = std::make_shared<penalized_contact_nonmatching_meshes_brick>
2400 (region1, region2,
false , option);
2402 tl.push_back(model::term_description(varname_u1, varname_u1,
true));
2403 tl.push_back(model::term_description(varname_u2, varname_u2,
true));
2404 tl.push_back(model::term_description(varname_u1, varname_u2,
true));
2405 tl.push_back(model::term_description(varname_u2, varname_u1,
true));
2407 model::varnamelist dl(1, dataname_r);
2410 case 2:
case 3: dl.push_back(dataname_lambda);
break;
2411 default: GMM_ASSERT1(
false,
"Penalized contact brick : invalid option");
2413 dl.push_back(dataname_friction_coeffs);
2414 if (dataname_alpha.size() > 0) {
2415 dl.push_back(dataname_alpha);
2416 if (dataname_wt1.size() > 0) {
2417 dl.push_back(dataname_wt1);
2418 if (dataname_wt2.size() > 0)
2419 dl.push_back(dataname_wt2);
2423 model::varnamelist vl(1, varname_u1);
2424 vl.push_back(varname_u2);
2426 return md.
add_brick(pbr, vl, dl, tl, model::mimlist(1, &mim), region1);
2431 void compute_integral_contact_area_and_force
2432 (model &md,
size_type indbrick, scalar_type &area,
2433 model_real_plain_vector &F) {
2435 pbrick pbr = md.brick_pointer(indbrick);
2436 const model::mimlist &ml = md.mimlist_of_brick(indbrick);
2437 const model::varnamelist &vl = md.varnamelist_of_brick(indbrick);
2438 const model::varnamelist &dl = md.datanamelist_of_brick(indbrick);
2439 size_type reg = md.region_of_brick(indbrick);
2441 GMM_ASSERT1(ml.size() == 1,
"Wrong size");
2443 if (pbr->brick_name() ==
"Integral contact with rigid obstacle brick" ||
2444 pbr->brick_name() ==
"Integral contact and friction with rigid obstacle brick") {
2445 integral_contact_rigid_obstacle_brick *p
2446 =
dynamic_cast<integral_contact_rigid_obstacle_brick *
>
2447 (
const_cast<virtual_brick *
>(pbr.get()));
2448 GMM_ASSERT1(p,
"Wrong type of brick");
2450 GMM_ASSERT1(vl.size() >= 2,
"Wrong size");
2451 const model_real_plain_vector &u = md.real_variable(vl[0]);
2452 const mesh_fem &mf_u = md.mesh_fem_of_variable(vl[0]);
2453 const model_real_plain_vector &lambda = md.real_variable(vl[1]);
2454 const mesh_fem &mf_lambda = md.mesh_fem_of_variable(vl[1]);
2455 GMM_ASSERT1(dl.size() >= 1,
"Wrong size");
2456 const model_real_plain_vector &obs = md.real_variable(dl[0]);
2457 const mesh_fem &mf_obs = md.mesh_fem_of_variable(dl[0]);
2460 area = asm_level_set_contact_area(*ml[0], mf_u, u, mf_obs, obs, reg, -1e-3,
2461 &mf_lambda, &lambda, 1e-1);
2463 gmm::resize(F, mf_u.nb_dof());
2465 (F, *ml[0], mf_u, mf_obs, obs, mf_lambda, lambda, reg);
2467 else if (pbr->brick_name() ==
"Integral penalized contact with rigid obstacle brick" ||
2468 pbr->brick_name() ==
"Integral penalized contact and friction with rigid "
2470 penalized_contact_rigid_obstacle_brick *p
2471 =
dynamic_cast<penalized_contact_rigid_obstacle_brick *
>
2472 (
const_cast<virtual_brick *
>(pbr.get()));
2473 GMM_ASSERT1(p,
"Wrong type of brick");
2474 GMM_ASSERT1(
false,
"Not implemented yet");
2476 else if (pbr->brick_name() ==
"Integral contact between nonmatching meshes brick" ||
2477 pbr->brick_name() ==
"Integral contact and friction between nonmatching "
2479 integral_contact_nonmatching_meshes_brick *p
2480 =
dynamic_cast<integral_contact_nonmatching_meshes_brick *
>
2481 (
const_cast<virtual_brick *
>(pbr.get()));
2482 GMM_ASSERT1(p,
"Wrong type of brick");
2484 GMM_ASSERT1(vl.size() == 3,
"Wrong size");
2485 const model_real_plain_vector &u1 = md.real_variable(vl[0]);
2486 const model_real_plain_vector &u2 = md.real_variable(vl[1]);
2487 const mesh_fem &mf_u1 = md.mesh_fem_of_variable(vl[0]);
2488 const mesh_fem &mf_u2 = md.mesh_fem_of_variable(vl[1]);
2489 const model_real_plain_vector &lambda = md.real_variable(vl[2]);
2490 const mesh_fem &mf_lambda = md.mesh_fem_of_variable(vl[2]);
2493 getfem::mesh_fem mf_u2_proj(mf_u1.linked_mesh(), mf_u1.linked_mesh().dim());
2494 mf_u2_proj.set_finite_element(mf_u1.linked_mesh().convex_index(), pfem_proj);
2496 std::vector<size_type> ind;
2497 mf_u2_proj.get_global_dof_index(ind);
2498 gmm::unsorted_sub_index SUBI(ind);
2501 size_type nbsub = mf_u2_proj.nb_basic_dof();
2502 model_real_plain_vector u2_proj(nbsub);
2504 if (mf_u2.is_reduced()) {
2505 gmm::csc_matrix<scalar_type> Esub(nbsub, nbdof2);
2506 gmm::copy(gmm::sub_matrix(mf_u2.extension_matrix(),
2507 SUBI, gmm::sub_interval(0, nbdof2)),
2512 gmm::copy(gmm::sub_vector(u2, SUBI), u2_proj);
2515 area = asm_nonmatching_meshes_contact_area
2516 (*ml[0], mf_u1, u1, mf_u2_proj, u2_proj, reg, -1e-3,
2517 &mf_lambda, &lambda, 1e-1);
2520 asm_nonmatching_meshes_normal_source_term
2521 (F, *ml[0], mf_u1, mf_u2_proj, mf_lambda, lambda, reg);
2525 else if (pbr->brick_name() ==
"Integral penalized contact between nonmatching meshes brick" ||
2526 pbr->brick_name() ==
"Integral penalized contact and friction between nonmatching "
2528 penalized_contact_nonmatching_meshes_brick *p
2529 =
dynamic_cast<penalized_contact_nonmatching_meshes_brick *
>
2530 (
const_cast<virtual_brick *
>(pbr.get()));
2531 GMM_ASSERT1(p,
"Wrong type of brick");
2532 GMM_ASSERT1(
false,
"Not implemented yet");
2544 (
model &md,
const mesh_im &mim,
const std::string &varname_u,
2545 const std::string &Neumannterm,
2546 const std::string &obs,
const std::string &dataname_gamma0,
2548 std::string dataname_friction_coeff,
2549 const std::string &dataname_alpha,
2550 const std::string &dataname_wt,
2553 std::string theta = std::to_string(theta_);
2554 ga_workspace workspace(md, ga_workspace::inherit::ALL);
2555 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
2556 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
2560 std::string gamma =
"(("+dataname_gamma0+
")/element_size)";
2561 std::string thetagamma =
"("+theta+
"/"+gamma+
")";
2562 std::string contact_normal =
"(-Normalized(Grad_"+obs+
"))";
2563 std::string gap =
"("+obs+
"-"+varname_u+
"."+contact_normal+
")";
2564 std::string Vs =
"("+varname_u +
2565 (dataname_wt.size() ?
"-("+dataname_wt+
"))" :
")");
2566 if (dataname_alpha.size()) Vs =
"(("+dataname_alpha+
")*"+Vs;
2567 if (dataname_friction_coeff.size() == 0)
2568 dataname_friction_coeff = std::string(
"0");
2570 std::string expr =
"Coulomb_friction_coupled_projection("+Neumannterm
2571 +
", "+contact_normal+
", "+Vs+
", "+gap+
", "+dataname_friction_coeff
2574 if (theta_ != scalar_type(0)) {
2575 std::string derivative_Neumann=workspace.extract_order1_term(varname_u);
2576 expr =
"-"+thetagamma+
"*("+Neumannterm+
").("+derivative_Neumann
2577 +
") + "+expr+thetagamma+
"*("+derivative_Neumann+
")";
2579 expr +=
"-Test_"+varname_u+
")";
2582 (md, mim, expr, region,
false,
false,
2583 "Nitsche contact with rigid obstacle");
2586 #ifdef EXPERIMENTAL_PURPOSE_ONLY
2596 size_type add_Nitsche_contact_with_rigid_obstacle_brick_modified_midpoint
2597 (model &md,
const mesh_im &mim,
const std::string &varname_u,
2598 const std::string &Neumannterm,
2599 const std::string &Neumannterm_wt,
2600 const std::string &obs,
const std::string &dataname_gamma0,
2602 std::string dataname_friction_coeff,
2603 const std::string &,
2604 const std::string &dataname_wt,
2607 std::string theta = std::to_string(theta_);
2608 ga_workspace workspace(md, ga_workspace::inherit::ALL);
2609 size_type order = workspace.add_expression(Neumannterm, mim, region, 1);
2610 GMM_ASSERT1(order == 0,
"Wrong expression of the Neumann term");
2612 std::string gamma =
"((1/("+dataname_gamma0+
"))*element_size)";
2613 std::string thetagamma =
"("+theta+
"*"+gamma+
")";
2614 std::string contact_normal =
"(-Normalized(Grad_"+obs+
"))";
2615 std::string gap =
"("+obs+
"-"+varname_u+
"."+contact_normal+
")";
2616 std::string gap_wt =
"("+obs+
"-"+dataname_wt+
"."+contact_normal+
")";
2617 std::string Vs =
"("+varname_u +
2618 (dataname_wt.size() ?
"-("+dataname_wt+
"))" :
")");
2619 std::string Vs_wt = Vs;
2620 if (dataname_friction_coeff.size() == 0)
2621 dataname_friction_coeff = std::string(
"0");
2623 std::string Pg_wt =
"(-"+gamma+
"*("+Neumannterm_wt+
")."+contact_normal+
2627 =
"((Heaviside("+Pg_wt+
")*Coulomb_friction_coupled_projection("
2628 +Neumannterm+
", "+contact_normal+
", "+Vs+
", "+gap+
", "
2629 +dataname_friction_coeff+
",1/"+gamma+
"))+"
2630 +
"((1-Heaviside("+Pg_wt+
"))*(Coulomb_friction_coupled_projection("
2631 +Neumannterm_wt+
", "+contact_normal+
", "+Vs_wt+
", "+gap_wt+
", "
2632 +dataname_friction_coeff+
",1/"+gamma+
")*0.5+"
2633 "Coulomb_friction_coupled_projection(2*("
2634 +Neumannterm+
")-("+Neumannterm_wt+
"), "+contact_normal
2635 +
", 2*"+Vs+
"-"+Vs_wt+
", 2*"+gap+
"-"+gap_wt+
", "
2636 +dataname_friction_coeff+
",1/"+gamma+
")*0.5))).(";
2638 if (theta_ != scalar_type(0)) {
2639 std::string derivative_Neumann=workspace.extract_order1_term(varname_u);
2640 expr =
"-"+thetagamma+
"*("+Neumannterm+
").("+derivative_Neumann
2641 +
") + "+expr+thetagamma+
"*("+derivative_Neumann+
")";
2643 expr +=
"-Test_"+varname_u+
")";
2646 (md, mim, expr, region,
false,
false,
2647 "Nitsche contact with rigid obstacle");
region-tree for window/point search on a set of rectangles.
Generic assembly of vectors, matrices.
void push_vec(VEC &v)
Add a new output vector.
void push_mat(const MAT &m)
Add a new output matrix (fake const version..)
void assembly(const mesh_region ®ion=mesh_region::all_convexes())
do the assembly on the specified region (boundary or set of convexes)
void push_nonlinear_term(pnonlinear_elem_term net)
Add a new non-linear term.
void push_mi(const mesh_im &im_)
Add a new mesh_im.
void push_mf(const mesh_fem &mf_)
Add a new mesh_fem.
Describe a finite element method linked to a mesh.
Describe an integration method linked to a mesh.
`‘Model’' variables store the variables, the data and the description of a model.
size_type add_brick(pbrick pbr, const varnamelist &varnames, const varnamelist &datanames, const termlist &terms, const mimlist &mims, size_type region)
Add a brick to the model.
FEM which projects a mesh_fem on a different mesh.
void copy(const L1 &l1, L2 &l2)
*/
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
void clear(L &l)
clear (fill with zeros) a vector or matrix.
void resize(V &v, size_type n)
*/
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*/
strongest_value_type< V1, V2 >::value_type vect_sp(const V1 &v1, const V2 &v2)
*/
computation of the condition number of dense matrices.
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
size_t size_type
used as the common size type in the library
size_type alpha(short_type n, short_type d)
Return the value of which is the number of monomials of a polynomial of variables and degree .
GEneric Tool for Finite Element Methods.
size_type add_integral_contact_between_nonmatching_meshes_brick(model &md, const mesh_im &mim, const std::string &varname_u1, const std::string &varname_u2, const std::string &multname_n, const std::string &dataname_r, size_type region1, size_type region2, int option=1)
Add a frictionless contact condition between nonmatching meshes to the model.
size_type APIDECL add_nonlinear_term(model &md, const mesh_im &mim, const std::string &expr, size_type region=size_type(-1), bool is_sym=false, bool is_coercive=false, const std::string &brickname="")
Add a nonlinear term given by the weak form language expression expr which will be assembled in regio...
size_type add_integral_contact_with_rigid_obstacle_brick(model &md, const mesh_im &mim, const std::string &varname_u, const std::string &multname_n, const std::string &dataname_obs, const std::string &dataname_r, size_type region, int option=1)
Add a frictionless contact condition with a rigid obstacle to the model, which is defined in an integ...
size_type add_penalized_contact_with_rigid_obstacle_brick(model &md, const mesh_im &mim, const std::string &varname_u, const std::string &dataname_obs, const std::string &dataname_r, size_type region, int option=1, const std::string &dataname_lambda_n="")
Add a penalized contact frictionless condition with a rigid obstacle to the model.
void asm_level_set_normal_source_term(VEC &R, const mesh_im &mim, const getfem::mesh_fem &mf_u, const getfem::mesh_fem &mf_obs, const VEC &obs, const getfem::mesh_fem &mf_lambda, const VEC &lambda, const mesh_region &rg)
Specific assembly procedure for the use of an Uzawa algorithm to solve contact problems.
void del_projected_fem(pfem pf)
release a projected fem
size_type add_penalized_contact_between_nonmatching_meshes_brick(model &md, const mesh_im &mim, const std::string &varname_u1, const std::string &varname_u2, const std::string &dataname_r, size_type region1, size_type region2, int option=1, const std::string &dataname_lambda_n="")
Add a penalized contact frictionless condition between nonmatching meshes to the model.
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
std::shared_ptr< const virtual_brick > pbrick
type of pointer on a brick
pfem new_projected_fem(const mesh_fem &mf_source, const mesh_im &mim_target, size_type rg_source_=size_type(-1), size_type rg_target_=size_type(-1), dal::bit_vector blocked_dofs=dal::bit_vector(), bool store_val=true)
create a new projected FEM.
size_type add_Nitsche_contact_with_rigid_obstacle_brick(model &md, const mesh_im &mim, const std::string &varname_u, const std::string &Neumannterm, const std::string &expr_obs, const std::string &dataname_gamma0, scalar_type theta_, std::string dataexpr_friction_coeff, const std::string &dataname_alpha, const std::string &dataname_wt, size_type region)
Adds a contact condition with or without Coulomb friction on the variable varname_u and the mesh boun...