32 const std::string &Theta,
33 const std::string ¶m_E,
34 const std::string ¶m_nu,
35 const std::string ¶m_epsilon,
36 const std::string ¶m_kappa,
41 std::string test_U =
"Test_" + sup_previous_and_dot_to_varname(U);
42 std::string test_Theta =
"Test_" + sup_previous_and_dot_to_varname(Theta);
43 std::string proj_Theta = (variant == 2) ?
44 (
"Elementary_transformation("+Theta+
",_2D_rotated_RT0_projection__434)")
46 std::string proj_test_Theta = (variant == 2) ?
47 (
"Elementary_transformation("+test_Theta
48 +
",_2D_rotated_RT0_projection__434)") : test_Theta;
50 std::string D =
"(("+param_E+
")*pow("+param_epsilon+
51 ",3))/(12*(1-sqr("+param_nu+
")))";
52 std::string G =
"(("+param_E+
")*("+param_epsilon+
"))*("+param_kappa+
53 ")/(2*(1+("+param_nu+
")))";
54 std::string E_theta =
"(Grad_" + Theta +
"+(Grad_" + Theta +
")')/2";
55 std::string E_test_Theta=
"(Grad_"+test_Theta+
"+(Grad_"+test_Theta+
")')/2";
57 std::string expr_left =
58 D+
"*(( 1-("+param_nu+
"))*("+E_theta+
"):("+E_test_Theta+
")+("+param_nu+
59 ")*Trace("+E_theta+
")*Trace("+E_test_Theta+
"))";
61 std::string expr_right =
62 "("+G+
")*(Grad_"+U+
"-"+proj_Theta+
").Grad_"+test_U+
63 "-("+G+
")*(Grad_"+U+
"-"+proj_Theta+
")."+proj_test_Theta;
67 return add_linear_generic_assembly_brick
68 (md, mim, expr_left+
"+"+expr_right, region,
false,
false,
69 "Reissner-Mindlin plate model brick");
71 add_linear_generic_assembly_brick
72 (md, mim, expr_left, region,
false,
false,
73 "Reissner-Mindlin plate model brick, rotation term");
74 return add_linear_generic_assembly_brick
75 (md, mim_red, expr_right, region,
false,
false,
76 "Reissner-Mindlin plate model brick, transverse shear term");
79 return add_linear_generic_assembly_brick
80 (md, mim, expr_left+
"+"+expr_right, region,
false,
false,
81 "Reissner-Mindlin plate model brick");
83 default: GMM_ASSERT1(
false,
"Invalid variant for Reissner-Mindlin brick.");
93 const std::string &Ua,
94 const std::string &Theta,
95 const std::string &U3,
96 const std::string &Theta3,
97 const std::string ¶m_E,
98 const std::string ¶m_nu,
99 const std::string ¶m_epsilon,
103 std::string test_Ua =
"Test_" + sup_previous_and_dot_to_varname(Ua);
104 std::string test_U3 =
"Test_" + sup_previous_and_dot_to_varname(U3);
105 std::string test_Theta =
"Test_" + sup_previous_and_dot_to_varname(Theta);
106 std::string proj_Theta = (variant >= 2) ?
107 (
"Elementary_transformation("+Theta+
",_2D_rotated_RT0_projection__434)")
109 std::string proj_test_Theta = (variant >= 2) ?
110 (
"Elementary_transformation("+test_Theta+
",_2D_rotated_RT0_projection__434)")
112 std::string test_Theta3 =
"Test_" + sup_previous_and_dot_to_varname(Theta3);
113 std::string proj_Theta3 = (variant == 3) ?
114 (
"Elementary_transformation("+Theta3+
",_P0_projection__434)")
116 std::string proj_test_Theta3 = (variant == 3) ?
117 (
"Elementary_transformation("+test_Theta3+
",_P0_projection__434)")
119 std::string D1 =
"(("+param_E+
")*pow("+param_epsilon+
",3))/(12*(1+("+param_nu+
")))";
120 std::string D2 = D1+
"*("+param_nu+
")/(1-2*("+param_nu+
"))";
121 std::string D3 = D1+
"/2";
122 std::string G1 =
"(("+param_E+
")*("+param_epsilon+
"))/(1+("+param_nu+
"))";
123 std::string G2 = G1+
"*("+param_nu+
")/(1-2*("+param_nu+
"))";
124 std::string G3 = G1+
"/2";
126 std::string E_Ua =
"(Grad_" + Ua +
"+(Grad_" + Ua +
")')/2";
127 std::string E_test_Ua =
"(Grad_" + test_Ua +
"+(Grad_" + test_Ua +
")')/2";
128 std::string E_Theta =
"(Grad_" + Theta +
"+(Grad_" + Theta +
")')/2";
129 std::string E_test_Theta=
"(Grad_"+test_Theta+
"+(Grad_"+test_Theta+
")')/2";
131 std::string expr_no_coupled_1 = G1+
"*("+E_Ua+
"):("+E_test_Ua+
") + "+G1+
"*("+Theta3+
")*("+test_Theta3+
")";
132 std::string expr_no_coupled_2 = D1+
"*("+E_Theta+
"):("+E_test_Theta+
") + "+D2+
"*Trace(Grad_"+Theta+
")*Trace(Grad_"+test_Theta+
") + "+D3+
"*(Grad_"+Theta3+
").(Grad_"+test_Theta3+
")";
134 std::string expr_coupled_1 = G3+
"*(Grad_"+U3+
" + "+proj_Theta+
").(Grad_"+test_U3+
" + "+proj_test_Theta+
")";
135 std::string expr_coupled_2 = G2+
"*(Trace("+E_Ua+
") + "+proj_Theta3+
")*(Trace("+E_test_Ua+
") + "+proj_test_Theta3+
")";
139 add_nonlinear_generic_assembly_brick
140 (md, mim, expr_no_coupled_1+
"+"+expr_no_coupled_2, region,
false,
false,
141 "enriched Reissner-Mindlin plate model brick, no coupled");
142 return add_nonlinear_generic_assembly_brick
143 (md, mim, expr_coupled_1+
"+"+expr_coupled_2, region,
false,
false,
144 "enriched Reissner-Mindlin plate model brick, coupled");
146 add_nonlinear_generic_assembly_brick
147 (md, mim, expr_no_coupled_1+
"+"+expr_no_coupled_2, region,
false,
false,
148 "enriched Reissner-Mindlin plate model brick, no coupled");
149 add_nonlinear_generic_assembly_brick
150 (md, mim_red1, expr_coupled_1, region,
false,
false,
151 "enriched Reissner-Mindlin plate model brick, coupled MR");
152 return add_nonlinear_generic_assembly_brick
153 (md, mim_red2, expr_coupled_2, region,
false,
false,
154 "enriched Reissner-Mindlin plate model brick, coupled eMR");
157 add_nonlinear_generic_assembly_brick
158 (md, mim, expr_no_coupled_1+
"+"+expr_no_coupled_2, region,
false,
false,
159 "enriched Reissner-Mindlin plate model brick, no coupled");
160 add_nonlinear_generic_assembly_brick
161 (md, mim, expr_coupled_1, region,
false,
false,
162 "enriched Reissner-Mindlin plate model brick, coupled MR");
163 return add_nonlinear_generic_assembly_brick
164 (md, mim_red2, expr_coupled_2, region,
false,
false,
165 "enriched Reissner-Mindlin plate model brick, coupled eMR");
169 add_nonlinear_generic_assembly_brick
170 (md, mim, expr_no_coupled_1+
"+"+expr_no_coupled_2, region,
false,
false,
171 "enriched Reissner-Mindlin plate model brick, no coupled");
172 add_nonlinear_generic_assembly_brick
173 (md, mim, expr_coupled_1, region,
false,
false,
174 "enriched Reissner-Mindlin plate model brick, coupled MR");
175 return add_nonlinear_generic_assembly_brick
176 (md, mim, expr_coupled_2, region,
false,
false,
177 "enriched Reissner-Mindlin plate model brick, coupled eMR");
179 default: GMM_ASSERT1(
false,
" testInvalid variant for enriched Reissner-Mindlin brick.");
193 class _2D_Rotated_RT0_projection_transformation
194 :
public virtual_elementary_transformation {
198 virtual void give_transformation(
const mesh_fem &mf,
const mesh_fem &mf2,
201 THREAD_SAFE_STATIC base_matrix M_old;
202 THREAD_SAFE_STATIC
pfem pf_old =
nullptr;
204 GMM_ASSERT1(&mf == &mf2,
205 "This transformation works on identical fems only");
208 pfem pf1 = mf.fem_of_element(cv);
210 GMM_ASSERT1(pf1->dim() == 2,
"This projection is only defined "
211 "for two-dimensional elements");
214 bool simplex =
false;
217 }
else if (pf1->ref_convex(cv)
221 GMM_ASSERT1(
false,
"Cannot adapt the method for such an element.");
224 if (pf1 == pf_old && pf1->is_equivalent() && M.size() == M_old.size()) {
229 std::stringstream fem_desc;
230 fem_desc <<
"FEM_RT0" << (simplex ?
"":
"Q") <<
"(" << N <<
")";
234 size_type degree = pf1->estimated_degree() + pf2->estimated_degree();
236 papprox_integration pim
240 size_type ndof1 = pf1->nb_dof(cv) * qmult;
242 base_matrix M1(ndof1, ndof1), M2(ndof2, ndof2), B(ndof1, ndof2);
243 base_matrix aux0(ndof1, ndof1), aux1(ndof1, ndof2), aux2(ndof1, ndof2);
244 base_matrix aux3(ndof2, ndof2);
248 bgeot::vectors_to_base_matrix(G, mf.linked_mesh().points_of_convex(cv));
249 fem_interpolation_context ctx1(pgt, pf1, base_node(N), G, cv);
250 fem_interpolation_context ctx2(pgt, pf2, base_node(N), G, cv);
253 base_matrix tv1, tv2;
255 for (
size_type i = 0; i < pim->nb_points_on_convex(); ++i) {
257 scalar_type coeff = pim->coeff(i);
258 ctx1.set_xref(pim->point(i));
259 ctx2.set_xref(pim->point(i));
260 pf1->real_base_value(ctx1, t1);
261 vectorize_base_tensor(t1, tv1, ndof1, pf1->target_dim(), N);
262 pf2->real_base_value(ctx2, t2);
263 vectorize_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
264 for (
size_type j = 0; j < ndof2; ++j) std::swap(tv2(j,0), tv2(j,1));
266 gmm::mult(tv1, gmm::transposed(tv1), aux0);
267 gmm::add(gmm::scaled(aux0, coeff), M1);
268 gmm::mult(tv2, gmm::transposed(tv2), aux3);
269 gmm::add(gmm::scaled(aux3, coeff), M2);
270 gmm::mult(tv1, gmm::transposed(tv2), aux1);
271 gmm::add(gmm::scaled(aux1, coeff), B);
278 gmm::mult(M1, B, aux1);
279 gmm::mult(aux1, M2, aux2);
280 GMM_ASSERT1(gmm::mat_nrows(M) == ndof1,
281 "Element not convenient for projection");
282 gmm::mult(aux2, gmm::transposed(B), M);
283 gmm::clean(M, 1E-15);
284 M_old = M; pf_old = pf1;
289 pelementary_transformation
290 p = std::make_shared<_2D_Rotated_RT0_projection_transformation>();
297 class _P0_projection_transformation
298 :
public virtual_elementary_transformation {
302 virtual void give_transformation(
const mesh_fem &mf,
const mesh_fem &mf2,
305 THREAD_SAFE_STATIC base_matrix M_old;
306 THREAD_SAFE_STATIC
pfem pf_old =
nullptr;
308 GMM_ASSERT1(&mf == &mf2,
309 "This transformation works on identical fems only");
312 pfem pf1 = mf.fem_of_element(cv);
313 size_type N = mf.get_qdim(), d = pf1->dim();
318 bool simplex =
false;
321 }
else if (pf1->ref_convex(cv)
325 GMM_ASSERT1(
false,
"Cannot adapt the method for such an element.");
328 if (pf1 == pf_old && pf1->is_equivalent() && M.size() == M_old.size()) {
333 std::stringstream fem_desc;
334 fem_desc <<
"FEM_" << (simplex ?
"PK":
"QK") <<
"(" << d <<
"," << 0 <<
")";
338 size_type degree = pf1->estimated_degree() + pf2->estimated_degree();
340 papprox_integration pim
344 size_type ndof1 = pf1->nb_dof(cv) * qmult;
345 size_type ndof2 = pf2->nb_dof(0) * qmult;
346 base_matrix M1(ndof1, ndof1), M2(ndof2, ndof2), B(ndof1, ndof2);
347 base_matrix aux0(ndof1, ndof1), aux1(ndof1, ndof2), aux2(ndof1, ndof2);
348 base_matrix aux3(ndof2, ndof2);
352 bgeot::vectors_to_base_matrix(G, mf.linked_mesh().points_of_convex(cv));
353 fem_interpolation_context ctx1(pgt, pf1, base_node(d), G, cv);
354 fem_interpolation_context ctx2(pgt, pf2, base_node(d), G, cv);
357 base_matrix tv1, tv2;
359 for (
size_type i = 0; i < pim->nb_points_on_convex(); ++i) {
361 scalar_type coeff = pim->coeff(i);
362 ctx1.set_xref(pim->point(i));
363 ctx2.set_xref(pim->point(i));
364 pf1->real_base_value(ctx1, t1);
365 vectorize_base_tensor(t1, tv1, ndof1, pf1->target_dim(), N);
366 pf2->real_base_value(ctx2, t2);
367 vectorize_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
370 gmm::mult(tv1, gmm::transposed(tv1), aux0);
371 gmm::add(gmm::scaled(aux0, coeff), M1);
372 gmm::mult(tv2, gmm::transposed(tv2), aux3);
373 gmm::add(gmm::scaled(aux3, coeff), M2);
374 gmm::mult(tv1, gmm::transposed(tv2), aux1);
375 gmm::add(gmm::scaled(aux1, coeff), B);
382 gmm::mult(M1, B, aux1);
383 gmm::mult(aux1, M2, aux2);
384 GMM_ASSERT1(gmm::mat_nrows(M) == ndof1,
385 "Element not convenient for projection");
386 gmm::mult(aux2, gmm::transposed(B), M);
387 gmm::clean(M, 1E-15);
388 M_old = M; pf_old = pf1;
396 pelementary_transformation
397 p = std::make_shared<_P0_projection_transformation>();