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;
 
   68         (md, mim, expr_left+
"+"+expr_right, region, 
false, 
false,
 
   69          "Reissner-Mindlin plate model brick");
 
   72         (md, mim, expr_left, region, 
false, 
false,
 
   73          "Reissner-Mindlin plate model brick, rotation term");
 
   75         (md, mim_red, expr_right, region, 
false, 
false,
 
   76          "Reissner-Mindlin plate model brick, transverse shear term");
 
   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+
")";
 
  140         (md, mim, expr_no_coupled_1+
"+"+expr_no_coupled_2, region, 
false, 
false,
 
  141          "enriched Reissner-Mindlin plate model brick, no coupled");
 
  143         (md, mim, expr_coupled_1+
"+"+expr_coupled_2, region, 
false, 
false,
 
  144          "enriched Reissner-Mindlin plate model brick, coupled");
 
  147         (md, mim, expr_no_coupled_1+
"+"+expr_no_coupled_2, region, 
false, 
false,
 
  148          "enriched Reissner-Mindlin plate model brick, no coupled");
 
  150         (md, mim_red1, expr_coupled_1, region, 
false, 
false,
 
  151          "enriched Reissner-Mindlin plate model brick, coupled MR");
 
  153         (md, mim_red2, expr_coupled_2, region, 
false, 
false,
 
  154          "enriched Reissner-Mindlin plate model brick, coupled eMR");
 
  158         (md, mim, expr_no_coupled_1+
"+"+expr_no_coupled_2, region, 
false, 
false,
 
  159          "enriched Reissner-Mindlin plate model brick, no coupled");
 
  161         (md, mim, expr_coupled_1, region, 
false, 
false,
 
  162          "enriched Reissner-Mindlin plate model brick, coupled MR");
 
  164         (md, mim_red2, expr_coupled_2, region, 
false, 
false,
 
  165          "enriched Reissner-Mindlin plate model brick, coupled eMR");  
 
  170         (md, mim, expr_no_coupled_1+
"+"+expr_no_coupled_2, region, 
false, 
false,
 
  171          "enriched Reissner-Mindlin plate model brick, no coupled");
 
  173         (md, mim, expr_coupled_1, region, 
false, 
false,
 
  174          "enriched Reissner-Mindlin plate model brick, coupled MR");
 
  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);
 
  280       GMM_ASSERT1(gmm::mat_nrows(M) == ndof1,
 
  281                   "Element not convenient for projection");
 
  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);
 
  384       GMM_ASSERT1(gmm::mat_nrows(M) == ndof1,
 
  385                   "Element not convenient for projection");
 
  388       M_old = M; pf_old = pf1;
 
  396     pelementary_transformation
 
  397       p = std::make_shared<_P0_projection_transformation>();
 
Describe an integration method linked to a mesh.
`‘Model’' variables store the variables, the data and the description of a model.
void add_elementary_transformation(const std::string &name, pelementary_transformation ptrans)
Add an elementary transformation to the model to be used with the generic assembly.
Reissner-Mindlin plate model brick.
void copy(const L1 &l1, L2 &l2)
*/
void clean(L &l, double threshold)
Clean a vector or matrix (replace near-zero entries with zeroes).
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*/
void add(const L1 &l1, L2 &l2)
*/
void lu_inverse(const DenseMatrixLU &LU, const Pvector &pvector, const DenseMatrix &AInv_)
Given an LU factored matrix, build the inverse of the matrix.
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
pfem fem_descriptor(const std::string &name)
get a fem descriptor from its string name.
pconvex_ref simplex_of_reference(dim_type nc, short_type K)
returns a simplex of reference of dimension nc and degree k
pconvex_ref parallelepiped_of_reference(dim_type nc, dim_type k)
parallelepiped of reference of dimension nc (and degree 1)
size_t size_type
used as the common size type in the library
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
GEneric Tool for Finite Element Methods.
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...
pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt, dim_type degree)
try to find an approximate integration method for the geometric transformation pgt which is able to i...
void add_P0_projection(model &md, std::string name)
Add the elementary transformation corresponding to the projection on P0 element.
size_type APIDECL add_linear_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="", bool return_if_nonlin=false)
Add a term given by the weak form language expression expr which will be assembled in region region a...
size_type add_Mindlin_Reissner_plate_brick(model &md, const mesh_im &mim, const mesh_im &mim_reduced, const std::string &u3, const std::string &Theta, const std::string ¶m_E, const std::string ¶m_nu, const std::string ¶m_epsilon, const std::string ¶m_kappa, size_type variant=size_type(2), size_type region=size_type(-1))
Add a term corresponding to the classical Reissner-Mindlin plate model for which u3 is the transverse...
void add_2D_rotated_RT0_projection(model &md, std::string name)
Add the elementary transformation corresponding to the projection on rotated RT0 element for two-dime...
size_type add_enriched_Mindlin_Reissner_plate_brick(model &md, const mesh_im &mim, const mesh_im &mim_reduced1, const mesh_im &mim_reduced2, const std::string &ua, const std::string &Theta, const std::string &u3, const std::string &Theta3, const std::string ¶m_E, const std::string ¶m_nu, const std::string ¶m_epsilon, size_type variant=size_type(3), size_type region=size_type(-1))
Add a term corresponding to the enriched Reissner-Mindlin plate model for which varname_ua is the mem...