[Getfem-commits] r4607 - in /trunk/getfem: interface/src/ interface/tests/matlab/ src/ src/getfem/

[Yves . Renard]

Author: renard
Date: Sun Apr 13 14:29:47 2014
New Revision: 4607

URL: http://svn.gna.org/viewcvs/getfem?rev=4607&view=rev
Log:
new nonlinear elsticity brick, in construction

Modified:
    trunk/getfem/interface/src/gf_model_set.cc
    trunk/getfem/interface/tests/matlab/demo_nonlinear_elasticity.m
    trunk/getfem/src/getfem/getfem_nonlinear_elasticity.h
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/src/getfem_nonlinear_elasticity.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=4607&r1=4606&r2=4607&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_model_set.cc  (original)
+++ trunk/getfem/interface/src/gf_model_set.cc  Sun Apr 13 14:29:47 2014
@@ -1464,8 +1464,9 @@
     depends on the law. It could be a short vector of constant values or a
     vector field described on a finite element method for variable
     coefficients. `region` is an optional mesh region on which the term
-    is added. If it is not specified, it is added on the whole mesh. Return the
-    brick index in the address@hidden/
+    is added. If it is not specified, it is added on the whole mesh.
+    This brick use the low-level generic assembly.
+    Returns the brick index in the address@hidden/
     sub_command
       ("add nonlinear elasticity brick", 4, 5, 0, 1,
        getfemint_mesh_im *gfi_mim = in.pop().to_getfemint_mesh_im();
@@ -1483,6 +1484,44 @@
        workspace().set_dependance(md, gfi_mim);
        out.pop().from_integer(int(ind));
        );
+
+    /address@hidden ind = ('add finite strain elasticity brick', @tmim mim, 
@str varname, @str constitutive_law, @str params[, @int region])
+    Add a nonlinear elasticity term to the model relatively to the
+    variable `varname`. `lawname` is the constitutive law which
+    could be 'SaintVenant Kirchhoff', 'Mooney Rivlin', 'Neo Hookean',
+    'Ciarlet Geymonat' or 'Generalized Blatz Ko'.
+    'Mooney Rivlin' and 'Neo Hookean' law names have to be preceeded with
+    the word 'Compressible' or 'Incompressible' to force using the
+    corresponding version.
+    The compressible version of these laws requires one additional material
+    coefficient.
+
+    IMPORTANT : if the variable is defined on a 2D mesh, the plane strain
+    approximation is automatically used.
+    `params` is a vector of parameters for the constitutive law. Its length
+    depends on the law. It could be a short vector of constant values or a
+    vector field described on a finite element method for variable
+    coefficients. `region` is an optional mesh region on which the term
+    is added. If it is not specified, it is added on the whole mesh.
+    This brick use the high-level generic assembly.
+    Returns the brick index in the address@hidden/
+    sub_command
+      ("add finite strain elasticity brick", 4, 5, 0, 1,
+       getfemint_mesh_im *gfi_mim = in.pop().to_getfemint_mesh_im();
+       size_type N = gfi_mim->mesh_im().linked_mesh().dim();
+       std::string varname = in.pop().to_string();
+       std::string lawname = in.pop().to_string();
+       std::string params = in.pop().to_string();
+       size_type region = size_type(-1);
+       if (in.remaining()) region = in.pop().to_integer();
+       size_type ind = config::base_index() +
+       add_finite_strain_elasticity_brick
+       (md->model(), gfi_mim->mesh_im(), varname, lawname, params, region);
+       workspace().set_dependance(md, gfi_mim);
+       out.pop().from_integer(int(ind));
+       );
+
+
 
     /address@hidden ind = ('add elastoplasticity brick', @tmim mim ,@str 
projname, @str varname, @str datalambda, @str datamu, @str datathreshold, @str 
datasigma[, @int region])
       Add a nonlinear elastoplastic term to the model relatively to the

Modified: trunk/getfem/interface/tests/matlab/demo_nonlinear_elasticity.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_nonlinear_elasticity.m?rev=4607&r1=4606&r2=4607&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_nonlinear_elasticity.m     
(original)
+++ trunk/getfem/interface/tests/matlab/demo_nonlinear_elasticity.m     Sun Apr 
13 14:29:47 2014
@@ -34,7 +34,7 @@
 lawname = 'SaintVenant Kirchhoff';
 params = [1;1];
 if (incompressible)
-    lawname = 'Mooney Rivlin';
+    lawname = 'Incompressible Mooney Rivlin';
     params = [1;1];
 end
 
@@ -85,7 +85,8 @@
 md=gf_model('real');
 gf_model_set(md, 'add fem variable', 'u', mfu);
 gf_model_set(md,'add initialized data','params', params);
-gf_model_set(md, 'add nonlinear elasticity brick', mim, 'u', lawname, 
'params');
+% gf_model_set(md, 'add nonlinear elasticity brick', mim, 'u', lawname, 
'params');
+gf_model_set(md, 'add finite strain elasticity brick', mim, 'u', lawname, 
'params');
 % gf_model_set(md, 'add nonlinear generic assembly brick', mim, ...
 %             'sqr(Trace(Green_Lagrangian(Id(meshdim)+Grad_u)))/8 + 
Norm_sqr(Green_Lagrangian(Id(meshdim)+Grad_u))/4');
 % gf_model_set(md, 'add nonlinear generic assembly brick', mim, ...

Modified: trunk/getfem/src/getfem/getfem_nonlinear_elasticity.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_nonlinear_elasticity.h?rev=4607&r1=4606&r2=4607&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_nonlinear_elasticity.h       (original)
+++ trunk/getfem/src/getfem/getfem_nonlinear_elasticity.h       Sun Apr 13 
14:29:47 2014
@@ -624,6 +624,18 @@
     gmm::copy(VMM, VM);
   }
   
+
+  /** Add a finite strain elasticity brick
+      to the model with respect to the variable
+      `varname` (the displacement).
+      For 2D meshes, switch automatically to plane strain elasticity.
+      High-level generic assembly version.
+  */
+  size_type add_finite_strain_elasticity_brick
+  (model &md, const mesh_im &mim, const std::string &varname,
+   std::string lawname, const std::string &params,
+   size_type region = size_type(-1));
+
 
   /** Add a nonlinear incompressibility term (for large strain elasticity)
       to the model with respect to the variable

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4607&r1=4606&r2=4607&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Sun Apr 13 14:29:47 2014
@@ -1348,7 +1348,7 @@
   static scalar_type ga_half_sqr_neg_part(scalar_type t)
   { return (t >= 0.) ? 0. : 0.5*t*t; }
   static scalar_type ga_sinc(scalar_type t) {// cardinal sine function sin(t)/t
-    if (gmm::abs(t) < 1E-5) {
+    if (gmm::abs(t) < 1E-4) {
       scalar_type t2 = t*t;
       return 1-t2/6.+ t2*t2/120.;
     } else {
@@ -1365,7 +1365,7 @@
   
   // Derivatives of predefined functions
   static scalar_type ga_der_sinc(scalar_type t) {
-    if (gmm::abs(t) < 1E-5) {
+    if (gmm::abs(t) < 1E-4) {
       scalar_type t2 = t*t;
       return  -t/3. + t*t2/30. -t*t2*t2/840.;
     } else {
@@ -1373,7 +1373,7 @@
     }
   }
   static scalar_type ga_der2_sinc(scalar_type t) {
-    if (gmm::abs(t) < 1E-5) {
+    if (gmm::abs(t) < 1E-4) {
       scalar_type t2 = t*t;
       return  -1./3. + t2/10. -t2*t2/168.;
     } else {

Modified: trunk/getfem/src/getfem_nonlinear_elasticity.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_nonlinear_elasticity.cc?rev=4607&r1=4606&r2=4607&view=diff
==============================================================================
--- trunk/getfem/src/getfem_nonlinear_elasticity.cc     (original)
+++ trunk/getfem/src/getfem_nonlinear_elasticity.cc     Sun Apr 13 14:29:47 2014
@@ -2071,14 +2071,59 @@
     return true;
   }
 
-
-
-
-
   static bool predef_operators_initialized = init_predef_operators();
 
 
 
+  size_type add_finite_strain_elasticity_brick
+  (model &md, const mesh_im &mim, const std::string &varname,
+   std::string lawname, const std::string &params,
+   size_type region) {
+    size_type N = mim.linked_mesh().dim();
+    GMM_ASSERT1(N >= 2 && N <= 3,
+                "Finite strain elasticity brick works only in 2D or 3D");
+
+    const mesh_fem *mf = md.pmesh_fem_of_variable(varname);
+    GMM_ASSERT1(mf, "Finite strain elasticity brick can only be applied on "
+                "fem variables");
+    size_type Q = mf->get_qdim();
+    GMM_ASSERT1(Q == N, "Finite strain elasticity brick can only be applied "
+                "on a fem variable having the same dimension than the mesh");
+
+    for (size_type i = 0; i < lawname.size(); ++i)
+      if (lawname[i] == ' ') lawname[i] = '_';
+
+    if (lawname.compare("SaintVenant_Kirchhoff") == 0) {
+      lawname = "Saint_Venant_Kirchhoff";
+    } else if (lawname.compare("Saint_Venant_Kirchhoff") == 0) {
+
+    } else if (lawname.compare("Generalized_Blatz_Ko") == 0) {
+      if (N == 2) lawname = "plane_strain_" + lawname;
+    } else if (lawname.compare("Ciarlet_Geymonat") == 0) {
+      if (N == 2) lawname = "plane_strain_" + lawname;
+    } else if (lawname.compare("Incompressible_Mooney_Rivlin") == 0) {
+      if (N == 2) lawname = "plane_strain_" + lawname;
+    } else if (lawname.compare("Compressible_Mooney_Rivlin") == 0) {
+      if (N == 2) lawname = "plane_strain_" + lawname;
+    } else if (lawname.compare("Incompressible_Neo_Hookean") == 0) {
+      if (N == 2) lawname = "plane_strain_" + lawname;
+    } else if (lawname.compare("Compressible_Neo_Hookean") == 0) {
+      if (N == 2) lawname = "plane_strain_" + lawname;
+    } else
+      GMM_ASSERT1(false, lawname << " is not a known hyperelastic law");
+
+    
+
+
+    std::string expr = "((Id(meshdim)+Grad_u)*(" + lawname
+      + "_sigma(Grad_"+varname+","+params+"))):Grad_Test_" + varname;
+
+    return add_nonlinear_generic_assembly_brick
+      (md, mim, expr, region, true, false,
+       "Finite strain elasticity brick for " + lawname + " law");
+  }
+
+
 
 }  /* end of namespace getfem.                                             */