[Getfem-commits] r4992 - in /trunk/getfem/interface/tests: matlab/ python/

[Yves . Renard]

Author: renard
Date: Wed May 13 17:40:54 2015
New Revision: 4992

URL: http://svn.gna.org/viewcvs/getfem?rev=4992&view=rev
Log:
Demo_Mindlin_Reissner in Python

Added:
    trunk/getfem/interface/tests/python/demo_Mindlin_Reissner_plate.py
Modified:
    trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m
    trunk/getfem/interface/tests/python/Makefile.am
    trunk/getfem/interface/tests/python/demo_crack.py

Modified: 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=4992&r1=4991&r2=4992&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m   
(original)
+++ trunk/getfem/interface/tests/matlab/demo_Mindlin_Reissner_plate.m   Wed May 
13 17:40:54 2015
@@ -1,4 +1,5 @@
-% Copyright (C) 20015-2015 FABRE Mathieu, SECK Mamadou, DALLERIT Valentin, 
Yves Renard
+% Copyright (C) 2015-2015 FABRE Mathieu, SECK Mamadou, DALLERIT Valentin,
+%                         Yves Renard
 %
 % This file is a part of GETFEM++
 %
@@ -24,12 +25,14 @@
 kappa     = 5/6;       % Shear correction factor
 f = -5*epsilon^3;      % Prescribed force on the top of the plate
 
-variant = 2;           % 0 : not reduced, 1 : with reduced integration, 2 : 
MITC reduction
+variant = 2;           % 0 : not reduced, 1 : with reduced integration,
+                       % 2 : MITC reduction
 quadrangles = true;    % Locking free only on quadrangle for the moment
 K = 1;                 % Degree of the finite element method
 with_Mindlin_brick = true; % Uses the Reissner-Mindlin predefined brick or not
-dirichlet_version = 1; % 0 = simplification, 1 = with multipliers, 2 = 
penalization
-
+dirichlet_version = 1; % 0 = simplification, 1 = with multipliers
+                       % 2 = penalization
+r = 1E8;               % Penalization parameter.
 plot_mesh = false;
 draw_solution = true;
 
@@ -42,11 +45,9 @@
   
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
+% Create a mesh_fem  for a 2D 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)));
@@ -59,15 +60,15 @@
   mim_reduced = gf_mesh_im(m, gf_integ('IM_TRIANGLE(1)'));
 end
 
-% detect the border of the mesh
+% Detect the border of the mesh and  assign it the boundary number 1
 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
 
+% Build the model
 md=gf_model('real');
 gf_model_set(md, 'add fem variable', 'u', mfu);
 gf_model_set(md, 'add fem variable', 'theta', mftheta);

Modified: trunk/getfem/interface/tests/python/Makefile.am
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/python/Makefile.am?rev=4992&r1=4991&r2=4992&view=diff
==============================================================================
--- trunk/getfem/interface/tests/python/Makefile.am     (original)
+++ trunk/getfem/interface/tests/python/Makefile.am     Wed May 13 17:40:54 2015
@@ -9,6 +9,7 @@
        demo_plasticity.py                              \
        demo_plate.py                                   \
        demo_static_contact.py                          \
+       demo_Mindlin_Reissner_plate.py                  \
        demo_step_by_step.py                            \
        demo_stokes_3D_tank.py                          \
        demo_stokes_3D_tank_draw.py                     \

Added: trunk/getfem/interface/tests/python/demo_Mindlin_Reissner_plate.py
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/python/demo_Mindlin_Reissner_plate.py?rev=4992&view=auto
==============================================================================
--- trunk/getfem/interface/tests/python/demo_Mindlin_Reissner_plate.py  (added)
+++ trunk/getfem/interface/tests/python/demo_Mindlin_Reissner_plate.py  Wed May 
13 17:40:54 2015
@@ -0,0 +1,172 @@
+#!/usr/bin/env python
+# -*- coding: UTF8 -*-
+# Python GetFEM++ interface
+#
+# Copyright (C) 2015-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 2.1 of the License,  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 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 demonstration.
+
+  This program is used to check that python-getfem is working. This is
+  also a good example of use of GetFEM++.
+
+"""
+import numpy as np
+import getfem as gf
+import os
+
+## Parameters
+Emodulus = 1           # Young Modulus
+nu       = 0.5         # Poisson Coefficient
+epsilon  = 0.001       # Plate thickness
+kappa     = 5/6        # Shear correction factor
+f = -5*pow(epsilon,3.) # Prescribed force on the top of the plate
+
+variant = 0            # 0 : not reduced, 1 : with reduced integration,
+                       # 2 : MITC reduction
+quadrangles = True     # Locking free only on quadrangle for the moment
+K = 1                  # Degree of the finite element method
+dirichlet_version = 1  # 0 = simplification, 1 = with multipliers,
+                       # 2 = penalization
+r = 1E8                # Penalization parameter.
+NX = 80                # Number of element per direction
+
+
+if (quadrangles):
+  m = gf.Mesh('cartesian', np.arange(0., 1., 1./NX),np.arange(0., 1., 1./NX))
+else:
+  m = gf.Mesh('import','structured',
+              'GT="GT_PK(2,1)";SIZES=[1,1];NOISED=0;NSUBDIV=[%d,%d];'
+              % (NX, NX))
+
+## Create a mesh_fem for a 2D vector field
+mftheta = gf.MeshFem(m, 2)
+mfu = gf.MeshFem(m, 1);
+mftheta.set_classical_fem(K)
+mfu.set_classical_fem(K)
+mim = gf.MeshIm(m, 6)
+mim_reduced = gf.MeshIm(m, 1)
+
+## Detect the border of the mesh and assign it the boundary number 1
+border = m.outer_faces()
+m.set_region(1, border)
+
+## Build the model
+md=gf.Model('real')
+md.add_fem_variable('u', mfu)
+md.add_fem_variable('theta', mftheta)
+md.add_initialized_data('E', Emodulus)
+md.add_initialized_data('nu', nu)
+md.add_initialized_data('epsilon', epsilon)
+md.add_initialized_data('kappa', kappa)
+
+md.add_Mindlin_Reissner_plate_brick(mim, mim_reduced, 'u', 'theta', 'E', 'nu', 
'epsilon', 'kappa', variant)
+
+md.add_initialized_data('VolumicData', f)
+md.add_source_term_brick(mim, 'u', 'VolumicData');
+md.add_initialized_data('DirichletData', 0);
+
+if (dirichlet_version == 0):
+    md.add_Dirichlet_condition_with_simplification('u', 1, 'DirichletData');   
+elif (dirichlet_version == 1):
+    md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu, 1, 
'DirichletData');
+elif (dirichlet_version == 2):
+    md.add_Dirichlet_condition_with_penalization(mim, 'u', r, 1, 
'DirichletData');
+
+print ('Number of Dof: %d' % md.nbdof()) 
+
+md.solve()
+U = md.variable('u');
+
+mfu.export_to_vtk('Deflection.vtk', U)
+print ('You can view solutions with for instance:\nmayavi2 -d Deflection.vtk 
-f WarpScalar -m Surface')
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+## Parameters
+PK=3; k = 1.;
+
+## Mesh and MeshFems
+m=Mesh('import','gid',filename);
+mfu=MeshFem(m,1);
+mfu.set_fem(Fem('FEM_PK(2,%d)' % (PK,)));
+mfd=MeshFem(m,1);
+mfd.set_fem(Fem('FEM_PK(2,%d)' % (PK,)));
+mim=MeshIm(m, Integ('IM_TRIANGLE(13)'));
+
+## Boundary selection
+P=m.pts(); # get list of mesh points coordinates
+Psqr=sum(P*P, 0);
+cobj=(Psqr < 1*1+1e-6);
+cout=(Psqr > 10*10-1e-2);
+pidobj=compress(cobj, range(0, m.nbpts()))
+pidout=compress(cout, range(0, m.nbpts()))
+fobj=m.faces_from_pid(pidobj)
+fout=m.faces_from_pid(pidout)
+ROBIN_BOUNDARY = 1
+DIRICHLET_BOUNDARY = 2
+m.set_region(DIRICHLET_BOUNDARY,fobj)
+m.set_region(ROBIN_BOUNDARY,fout)
+
+## Interpolate the exact solution on mfd (assuming it is a Lagrange fem)
+wave_expr = ('cos(%f*y+.2)+complex(0.,1.)*sin(%f*y+.2)' % (k,k));
+Uinc=mfd.eval(wave_expr,globals(),locals());
+
+## Model Bricks
+md = Model('complex')
+md.add_fem_variable('u', mfu)
+md.add_initialized_data('k', [k]);
+md.add_Helmholtz_brick(mim, 'u', 'k');
+md.add_initialized_data('Q', [complex(0.,1.)*k]);
+md.add_Fourier_Robin_brick(mim, 'u', 'Q', ROBIN_BOUNDARY);
+md.add_initialized_fem_data('DirichletData', mfd, Uinc);
+md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu, DIRICHLET_BOUNDARY,
+                                            'DirichletData');
+#md.add_Dirichlet_condition_with_penalization(mim, 'u', 1e12,
+#                                             DIRICHLET_BOUNDARY,
+#                                             'DirichletData');
+
+## Solving the problem
+md.solve();
+U = md.variable('u');
+
+if (not(make_check)):
+
+    sl=Slice(('none',), mfu, 8)
+    sl.export_to_vtk('wave.vtk', mfu, real(U), 'rWave',
+                     mfu, imag(U), 'iWave')
+
+    print 'You can view the solution with (for instance):'
+    print 'mayavi2 -d wave.vtk -f WarpScalar -m Surface'

Modified: trunk/getfem/interface/tests/python/demo_crack.py
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/python/demo_crack.py?rev=4992&r1=4991&r2=4992&view=diff
==============================================================================
--- trunk/getfem/interface/tests/python/demo_crack.py   (original)
+++ trunk/getfem/interface/tests/python/demo_crack.py   Wed May 13 17:40:54 2015
@@ -43,8 +43,8 @@
                                     #
 DIRICHLET  = 101                    #
                                     #
-Lambda = 1.25E10 # Lamé coefficient #
-Mu = 1.875E10    # Lamé coefficient #
+Lambda = 1.25E10 # Lame coefficient #
+Mu = 1.875E10    # Lame coefficient #
 #####################################
 
 # Global Functions: ###############################