[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] r4911 - in /trunk/getfem: interface/src/ interface/test
From: |
Yves . Renard |
Subject: |
[Getfem-commits] r4911 - in /trunk/getfem: interface/src/ interface/tests/matlab/ src/ src/getfem/ |
Date: |
Thu, 26 Mar 2015 13:16:45 -0000 |
Author: renard
Date: Thu Mar 26 14:16:45 2015
New Revision: 4911
URL: http://svn.gna.org/viewcvs/getfem?rev=4911&view=rev
Log:
Mindlin-Reissner plate model with MITC element
Added:
trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m
Modified:
trunk/getfem/interface/src/gf_model_set.cc
trunk/getfem/interface/tests/matlab/Makefile.am
trunk/getfem/src/getfem/getfem_linearized_plates.h
trunk/getfem/src/getfem_generic_assembly.cc
trunk/getfem/src/getfem_linearized_plates.cc
Modified: trunk/getfem/interface/src/gf_model_set.cc
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_model_set.cc?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_model_set.cc (original)
+++ trunk/getfem/interface/src/gf_model_set.cc Thu Mar 26 14:16:45 2015
@@ -1755,9 +1755,6 @@
out.pop().from_integer(int(ind));
);
-
-
-
/address@hidden ind = ('add Kirchhoff-Love plate brick', @tmim mim, @str
varname, @str dataname_D, @str dataname_nu [, @int region])
Add a bilaplacian brick on the variable
`varname` and on the mesh region `region`.
@@ -1927,6 +1924,55 @@
);
+ /address@hidden ind = ('add Reisner-Mindlin plate brick', @tmim mim, @tmim
mim_reduced, @str varname_U, @str varname_theta , @str param_E, @str param_nu,
@str param_epsilon, @str param_kappa [,@int variant [, @int region]])
+ Add a term corresponding to the classical Reissner-Mindlin plate
+ model for which `U` is the transverse displacement,
+ `Theta` the rotation of
+ fibers normal to the midplane, 'param_E' the Young Modulus,
+ `param_nu` the poisson ratio,
+ `param_epsilon` the plate thickness,
+ `param_kappa` the shear correction factor. Note that since this brick
+ uses the high level generic assembly language, the parameter can
+ be regular expression of this language.
+ There are three variants.
+ `variant = 0` corresponds to the an
+ unreduced formulation and in that case only the integration
+ method `mim` is used. Practically this variant is not usable since
+ it is subject to a strong locking phenomenon.
+ `variant = 1` corresponds to a reduced integration where `mim` is
+ used for the rotation term and `mim_reduced` for the transverse
+ shear term. `variant = 2` (default) corresponds to the projection onto
+ a rotated RT0 element of the transverse shear term. For the moment, this
+ is adapted to quadrilateral only (because it is not sufficient to
+ remove the locking phenomenon on triangle elements). Note also that if
+ you use high order elements, the projection on RT0 will reduce the order
+ of the approximation.
+ Returns the brick index in the model.
+ @*/
+ sub_command
+ ("add Mindlin Reissner plate brick", 7, 9, 0, 1,
+ getfemint_mesh_im *gfi_mim = in.pop().to_getfemint_mesh_im();
+ getfemint_mesh_im *gfi_mim_reduced = in.pop().to_getfemint_mesh_im();
+ std::string varname_U = in.pop().to_string();
+ std::string varname_theta = in.pop().to_string();
+ std::string param_E = in.pop().to_string();
+ std::string param_nu = in.pop().to_string();
+ std::string param_epsilon = in.pop().to_string();
+ std::string param_kapa = in.pop().to_string();
+ size_type variant = size_type(2);
+ if (in.remaining()) variant = in.pop().to_integer();
+ size_type region = size_type(-1);
+ if (in.remaining()) region = in.pop().to_integer();
+ size_type ind = add_Mindlin_Reissner_plate_brick
+ (md->model(), gfi_mim->mesh_im(), gfi_mim_reduced->mesh_im(),
+ varname_U, varname_theta, param_E, param_nu, param_epsilon,
+ param_kapa, variant, region);
+ workspace().set_dependance(md, gfi_mim);
+ out.pop().from_integer(int(ind));
+ );
+
+
+
/address@hidden ind = ('add mass brick', @tmim mim, @str varname[, @str
dataname_rho[, @int region]])
Add mass term to the model relatively to the variable `varname`.
If specified, the data `dataname_rho` should contain the
Modified: trunk/getfem/interface/tests/matlab/Makefile.am
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/Makefile.am?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/Makefile.am (original)
+++ trunk/getfem/interface/tests/matlab/Makefile.am Thu Mar 26 14:16:45 2015
@@ -25,6 +25,7 @@
check_plasticity.m \
demo_elasticity.m \
demo_bilaplacian.m \
+ demo_Mindlin_Reissner_plate.m \
demo_laplacian.m \
demo_periodic_laplacian.m \
demo_laplacian1D.m \
Added: trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m?rev=4911&view=auto
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m (added)
+++ trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m Thu Mar
26 14:16:45 2015
@@ -0,0 +1,114 @@
+% Copyright (C) 20015-2015 FABRE Mathieu, SECK Mamadou, DALLERIT Valentin,
Yves Renard
+%
+% This file is a part of GETFEM++
+%
+% Getfem++ is free software; you can redistribute it and/or modify it
+% under the terms of the GNU Lesser General Public License as published
+% by the Free Software Foundation; either version 3 of the License, or
+% (at your option) any later version along with the GCC Runtime Library
+% Exception either version 3.1 or (at your option) any later version.
+% This program is distributed in the hope that it will be useful, but
+% WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+% License and GCC Runtime Library Exception for more details.
+% You should have received a copy of the GNU Lesser General Public License
+% along with this program; if not, write to the Free Software Foundation,
+% Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+
+% Simple supported Mindlin-Reissner plate
+
+
+Emodulus = 1; % Young Modulus
+nu = 0.5; % Poisson Coefficient
+epsilon = 0.001; % Plate thickness
+kappa = 5/6; % Shear correction factor
+f = -10*epsilon^3; % Prescribed force on the top of the plate
+
+dirichlet_version = 1; % 0 = simplification, 1 = with multipliers, 2 =
penalization
+quadrangles = true; % Locking free only on quadrangle for the moment
+K = 1; % Degree of the finite element method
+variant = 1; % 0 : not reduced, 1 : with reduced integration, 2 :
MITC reduction
+with_Mindlin_brick = true; % Uses the Reissner-Mindlin predefined brick or not
+
+plot_mesh = false;
+draw_solution = true;
+
+% trace on;
+gf_workspace('clear all');
+NX = 20;
+if (quadrangles)
+ m = gf_mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
+else
+
m=gf_mesh('import','structured',sprintf('GT="GT_PK(2,1)";SIZES=[1,1];NOISED=0;NSUBDIV=[%d,%d];',
NX, NX));
+end
+
+% Create a mesh_fem of for a 2 dimension vector field
+mftheta = gf_mesh_fem(m,2);
+mfu = gf_mesh_fem(m,1);
+% Assign the QK or PK fem to all convexes of the mesh_fem, and define an
+% integration method
+if (quadrangles)
+ gf_mesh_fem_set(mftheta,'fem',gf_fem(sprintf('FEM_QK(2,%d)', K)));
+ gf_mesh_fem_set(mfu,'fem',gf_fem(sprintf('FEM_QK(2,%d)', K)));
+ mim = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,6)'));
+ mim_reduced = gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2,1)'));
+else
+ gf_mesh_fem_set(mftheta,'fem',gf_fem(sprintf('FEM_PK(2,%d)', K)));
+ gf_mesh_fem_set(mfu,'fem',gf_fem(sprintf('FEM_PK(2,%d)', K)));
+ mim = gf_mesh_im(m, gf_integ('IM_TRIANGLE(6)'));
+ mim_reduced = gf_mesh_im(m, gf_integ('IM_TRIANGLE(1)'));
+end
+
+% detect the border of the mesh
+border = gf_mesh_get(m,'outer faces');
+% mark it as boundary #1
+gf_mesh_set(m, 'boundary', 1, border);
+if (plot_mesh)
+ gf_plot_mesh(m, 'regions', [1]); % the boundary edges appears in red
+ pause(1);
+end
+
+md=gf_model('real');
+gf_model_set(md, 'add elementary rotated RT0 projection', 'RT0_projection');
+gf_model_set(md, 'add fem variable', 'u', mfu);
+gf_model_set(md, 'add fem variable', 'theta', mftheta);
+gf_model_set(md, 'add initialized data', 'E', Emodulus);
+gf_model_set(md, 'add initialized data', 'nu', nu);
+gf_model_set(md, 'add initialized data', 'epsilon', epsilon);
+gf_model_set(md, 'add initialized data', 'kappa', kappa);
+
+
+if (with_Mindlin_brick)
+ gf_model_set(md, 'add Mindlin Reissner plate brick', mim, mim_reduced, 'u',
'theta', 'E', 'nu', 'epsilon', 'kappa', variant);
+else
+ gf_model_set(md, 'add linear generic assembly brick', mim,
'(E*epsilon*epsilon*epsilon*(1-nu)/(48 * (1 - nu*nu))) *
((Grad_theta+Grad_theta''):(Grad_Test_theta+Grad_Test_theta''))');
+ gf_model_set(md, 'add linear generic assembly brick', mim,
'(E*epsilon*epsilon*epsilon*nu/(12 * (1 - nu*nu))) *
(Trace(Grad_theta)*Trace(Grad_Test_theta))');
+ if (variant == 0)
+ gf_model_set(md, 'add linear generic assembly brick', mim,
'(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Grad_Test_u) +
(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Test_theta)');
+ elseif (variant == 1)
+ gf_model_set(md, 'add linear generic assembly brick', mim_reduced,
'(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Grad_Test_u) +
(E*kappa*epsilon/(1 + nu)) * ((Grad_u + theta).Test_theta)');
+ else
+ gf_model_set(md, 'add linear generic assembly brick', mim,
'(E*kappa*epsilon/(1 + nu)) * ((Grad_u +
Elementary_transformation(theta,RT0_projection)).Grad_Test_u) +
(E*kappa*epsilon/(1 + nu)) * ((Grad_u + Elementary_transformation(theta,
RT0_projection)).(Elementary_transformation(Test_theta, RT0_projection)))');
+ end
+end
+
+gf_model_set(md, 'add initialized data', 'VolumicData', f);
+
+gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');
+gf_model_set(md, 'add initialized data', 'DirichletData', 0);
+switch (dirichlet_version)
+ case 0,
+ gf_model_set(md, 'add Dirichlet condition with simplification', 'u', 1,
'DirichletData');
+ case 1,
+ gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u',
mfu, 1, 'DirichletData');
+ case 2,
+ gf_model_set(md, 'add Dirichlet condition with penalization', mim, 'u', r,
1, 'DirichletData');
+end
+gf_model_get(md, 'solve');
+U = gf_model_get(md, 'variable', 'u');
+
+if (draw)
+ gf_plot(mfu,U,'mesh','off', 'zplot', 'on');
+ colorbar; title('computed solution');
+end
+
Modified: trunk/getfem/src/getfem/getfem_linearized_plates.h
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_linearized_plates.h?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_linearized_plates.h (original)
+++ trunk/getfem/src/getfem/getfem_linearized_plates.h Thu Mar 26 14:16:45 2015
@@ -49,6 +49,45 @@
*/
void add_2D_rotated_RT0_projection(model &md, std::string name);
+
+ /** Add a term corresponding to the classical Reissner-Mindlin plate
+ model for which `U` is the transverse displacement,
+ `Theta` the rotation of
+ fibers normal to the midplane, 'param_E' the Young Modulus,
+ `param_nu` the poisson ratio,
+ `param_epsilon` the plate thickness,
+ `param_kappa` the shear correction factor. Note that since this brick
+ uses the high level generic assembly language, the parameter can
+ be regular expression of this language.
+ There are three variants.
+ `variant = 0` corresponds to the an
+ unreduced formulation and in that case only the integration
+ method `mim` is used. Practically this variant is not usable since
+ it is subject to a strong locking phenomenon.
+ `variant = 1` corresponds to a reduced integration where `mim` is
+ used for the rotation term and `mim_reduced` for the transverse
+ shear term. `variant = 2` (default) corresponds to the projection onto
+ a rotated RT0 element of the transverse shear term. For the moment, this
+ is adapted to quadrilateral only (because it is not sufficient to
+ remove the locking phenomenon on triangle elements). Note also that if
+ you use high order elements, the projection on RT0 will reduce the order
+ of the approximation.
+ Returns the brick index in the model.
+ */
+ size_type add_Mindlin_Reissner_plate_brick
+ (model &md, const mesh_im &mim, const mesh_im &mim_reduced,
+ const std::string &U,
+ 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));
+
+
+
+
+
+
+
/* ******************************************************************** */
/* Linear plate specific assembly procedures. */
/* ******************************************************************** */
Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Thu Mar 26 14:16:45 2015
@@ -232,14 +232,17 @@
static void ga_throw_error_msg(const std::string &expr, size_type pos,
const std::string &msg) {
+ int lenght_before = 40, lenght_after = 40;
if (expr.size()) {
- int first = std::max(0, int(pos)-40);
- int last = std::min(int(pos)+20, int(expr.size()));
- if (last - first < 60)
- first = std::max(0, int(pos)-40-(60-last+first));
- if (last - first < 60)
- last = std::min(int(pos)+20+(60-last+first),int(expr.size()));
-
+ int first = std::max(0, int(pos)-lenght_before);
+ int last = std::min(int(pos)+lenght_after, int(expr.size()));
+ if (last - first < lenght_before+lenght_after)
+ first = std::max(0, int(pos)-lenght_before
+ -(lenght_before+lenght_after-last+first));
+ if (last - first < lenght_before+lenght_after)
+ last = std::min(int(pos)+lenght_after
+ +(lenght_before+lenght_after-last+first),
+ int(expr.size()));
if (first > 0) cerr << "...";
cerr << expr.substr(first, last-first);
if (last < int(expr.size())) cerr << "...";
@@ -6215,8 +6218,8 @@
{ test = 1; name = name.substr(5); }
if (!(workspace.variable_exists(name)))
- ga_throw_error(expr, pnode->pos,
- "Unknown variable, function, operator or data");
+ ga_throw_error(expr, pnode->pos, "Unknown variable, function, "
+ "operator or data " + name);
if (pnode->der1)
ga_throw_error(expr, pnode->pos, "Derivative is for functions or "
Modified: trunk/getfem/src/getfem_linearized_plates.cc
URL:
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_linearized_plates.cc?rev=4911&r1=4910&r2=4911&view=diff
==============================================================================
--- trunk/getfem/src/getfem_linearized_plates.cc (original)
+++ trunk/getfem/src/getfem_linearized_plates.cc Thu Mar 26 14:16:45 2015
@@ -1,7 +1,9 @@
/* -*- c++ -*- (enables emacs c++ mode) */
/*===========================================================================
- Copyright (C) 2004-2015 Yves Renard
+ Copyright (C) 2004-2015 Yves Renard, Jeremie Lasry, Mathieu Fabre
+ SECK Mamadou, DALLERIT Valentin
+
This file is a part of GETFEM++
@@ -26,6 +28,75 @@
namespace getfem {
+
+ size_type add_Mindlin_Reissner_plate_brick(model &md,
+ const mesh_im & mim,
+ const mesh_im & mim_red,
+ const std::string &U,
+ 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 region) {
+
+
+ std::string test_U = "Test_" + sup_previous_and_dot_to_varname(U);
+ std::string test_Theta = "Test_" + sup_previous_and_dot_to_varname(Theta);
+ std::string proj_Theta = (variant == 2) ?
+ ("Elementary_transformation("+Theta+",_2D_rotated_RT0_projection__434)")
+ : Theta;
+ std::string proj_test_Theta = (variant == 2) ?
+ ("Elementary_transformation("+test_Theta
+ +",_2D_rotated_RT0_projection__434)") : test_Theta;
+
+ std::string D = "(("+param_E+")*pow("+param_epsilon+
+ ",3))/(12*(1-sqr("+param_nu+")))";
+ std::string G = "(("+param_E+")*("+param_epsilon+"))*("+param_kappa+
+ ")/(2*(1+("+param_nu+")))";
+ std::string E_theta = "(Grad_" + Theta + "+(Grad_" + Theta + ")')/2";
+ std::string E_test_Theta="(Grad_"+test_Theta+"+(Grad_"+test_Theta+")')/2";
+
+ std::string expr_left =
+ D+"*(( 1-("+param_nu+"))*("+E_theta+"):("+E_test_Theta+")+("+param_nu+
+ ")*Trace("+E_theta+")*Trace("+E_test_Theta+"))";
+
+ std::string expr_right =
+ "("+G+")*(Grad_"+U+"-"+proj_Theta+").Grad_"+test_U+
+ "-("+G+")*(Grad_"+U+"-"+proj_Theta+")."+proj_test_Theta;
+
+ switch(variant) {
+ case 0: // Without reduction
+ return add_linear_generic_assembly_brick
+ (md, mim, expr_left+"+"+expr_right, region, false, false,
+ "Reissner-Mindlin plate model brick");
+ case 1: // With reduced integration
+ add_linear_generic_assembly_brick
+ (md, mim, expr_left, region, false, false,
+ "Reissner-Mindlin plate model brick, rotation term");
+ return add_linear_generic_assembly_brick
+ (md, mim_red, expr_right, region, false, false,
+ "Reissner-Mindlin plate model brick, transverse shear term");
+ case 2: // Variant with projection on rotated RT0
+ add_2D_rotated_RT0_projection(md, "_2D_rotated_RT0_projection__434");
+ return add_linear_generic_assembly_brick
+ (md, mim, expr_left+"+"+expr_right, region, false, false,
+ "Reissner-Mindlin plate model brick");
+ break;
+ default: GMM_ASSERT1(false, "Invalid variant for Reissner-Mindlin brick.");
+ }
+ return size_type(-1);
+ }
+
+
+
+
+
+
+
+
+
// For the moment, only projection onto rotated RT0 element in dimension 2
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Getfem-commits] r4911 - in /trunk/getfem: interface/src/ interface/tests/matlab/ src/ src/getfem/,
Yves . Renard <=