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

[Yves . Renard]

Author: renard
Date: Fri Jul  4 14:50:17 2014
New Revision: 4702

URL: http://svn.gna.org/viewcvs/getfem?rev=4702&view=rev
Log:
work in progress

Modified:
    trunk/getfem/interface/tests/matlab/demo_large_sliding_contact.m
    trunk/getfem/src/getfem/getfem_generic_assembly.h
    trunk/getfem/src/getfem_generic_assembly.cc

Modified: trunk/getfem/interface/tests/matlab/demo_large_sliding_contact.m
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/tests/matlab/demo_large_sliding_contact.m?rev=4702&r1=4701&r2=4702&view=diff
==============================================================================
--- trunk/getfem/interface/tests/matlab/demo_large_sliding_contact.m    
(original)
+++ trunk/getfem/interface/tests/matlab/demo_large_sliding_contact.m    Fri Jul 
 4 14:50:17 2014
@@ -29,14 +29,14 @@
 clambda2 = 1.; cmu2 = 1.;   % Elasticity parameters
 r = 0.1;                    % Augmentation parameter
 alpha = 1.0;                % Alpha coefficient for "sliding velocity"
-f_coeff = 0.0;              % Friction coefficient
-
-test_tangent_matrix = true;
+f_coeff = 1.0;              % Friction coefficient
+
+test_tangent_matrix = false;
 nonlinear_elasticity = false;
 max_iter = 50;
 draw_mesh = false;
 do_plot = true;
-generic_assembly_contact_brick = false;
+generic_assembly_contact_brick = true;
 
 switch(test_case)
   case {0,1}
@@ -250,11 +250,11 @@
   gf_model_set(md, 'add initialized data', 'r', r);
   gf_model_set(md, 'add initialized data', 'f', f_coeff);
   
-  % gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'-lambda1.Test_u1', CONTACT_BOUNDARY1); 
-  % gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'-(1/r)*lambda1.Test_lambda1', CONTACT_BOUNDARY1);
-  % gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'Interpolate_filter(contact_trans, 
(1/r)*Coulomb_friction_coupled_projection(lambda1, 
Transformed_unit_vector(Grad_u1, Normal), u1, 
(Interpolate(x,contact_trans)-x-u1).Transformed_unit_vector(Grad_u1, Normal), 
f, r).Test_lambda1, 2)', CONTACT_BOUNDARY1);
-  % gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'Interpolate_filter(contact_trans, 
(1/r)*Coulomb_friction_coupled_projection(lambda1, 
Transformed_unit_vector(Grad_u1, Normal), u1-Interpolate(u1,contact_trans), 
(x+u1-Interpolate(x,contact_trans)-Interpolate(u1,contact_trans)).Transformed_unit_vector(Grad_u1,
 Normal), f, r).Test_lambda1, 1)', CONTACT_BOUNDARY1);
-  gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'Interpolate_filter(contact_trans, 
(x-Interpolate(x,contact_trans)).Test_lambda1, 1)', CONTACT_BOUNDARY1);
+  gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'-lambda1.Test_u1', CONTACT_BOUNDARY1); 
+  gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'-lambda1.Test_lambda1', CONTACT_BOUNDARY1);
+  gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'Interpolate_filter(contact_trans, Coulomb_friction_coupled_projection(lambda1, 
Transformed_unit_vector(Grad_u1, Normal), u1, 
(Interpolate(x,contact_trans)-x-u1).Transformed_unit_vector(Grad_u1, Normal), 
f, r).Test_lambda1, 2)', CONTACT_BOUNDARY1);
+  gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'Interpolate_filter(contact_trans, Coulomb_friction_coupled_projection(lambda1, 
Transformed_unit_vector(Grad_u1, Normal), u1-Interpolate(u1,contact_trans), 
(Interpolate(x,contact_trans)+Interpolate(u1,contact_trans)-x-u1).Transformed_unit_vector(Grad_u1,
 Normal), f, r).Test_lambda1, 1)', CONTACT_BOUNDARY1);
+  % gf_model_set(md, 'add nonlinear generic assembly brick', mim1_contact, 
'Interpolate_filter(contact_trans, (Interpolate(x,contact_trans)).Test_lambda1, 
1)', CONTACT_BOUNDARY1);
     
 else
   mcff=gf_multi_contact_frame(md, N, release_dist, false, self_contact, 0.2, 
true, 0, false);
@@ -289,22 +289,22 @@
   gf_model_set(md, 'add integral large sliding contact brick raytrace', mcff, 
'r', 'f', 'alpha');
 end
 
-for nit=1:15
+for nit=1:150
   disp(sprintf('Iteration %d', nit));
 
   if (test_tangent_matrix) 
-    errmax = gf_model_get(md, 'test tangent matrix', 1E-8, 10, 0.00001);
+    errmax = gf_model_get(md, 'test tangent matrix', 1E-10, 10, 0.0001);
     % errmax = gf_model_get(md, 'test tangent matrix term', 'lambda1', 'u1', 
1E-8, 20, 0.00001);
     disp(sprintf('errmax = %g', errmax));
     if (errmax > 1E-3) error('bad tangent matrix'); end;
-    return;
-    pause(2);
-  end
-    
+    % return;
+    pause(0.1);
+  end
+
+  gf_model_get(md, 'solve', 'noisy', 'max_iter', max_iter, 'max_res', 
max_res); % , 'lsearch', 'simplest');
+
   
   if (do_plot)
-    gf_model_get(md, 'solve', 'noisy', 'max_iter', max_iter, 'max_res', 
max_res); % , 'lsearch', 'simplest');
-
     U1 = gf_model_get(md, 'variable', 'u1');
     if (nonlinear_elasticity)
       VM1 = gf_model_get(md, 'compute Von Mises or Tresca', 'u1', lawname, 
'params1', mfvm1);

Modified: trunk/getfem/src/getfem/getfem_generic_assembly.h
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem/getfem_generic_assembly.h?rev=4702&r1=4701&r2=4702&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_generic_assembly.h   (original)
+++ trunk/getfem/src/getfem/getfem_generic_assembly.h   Fri Jul  4 14:50:17 2014
@@ -149,7 +149,7 @@
       std::string interpolate_name_test1, interpolate_name_test2;
       const mesh_im *mim;
       const mesh *m;
-      mesh_region rg;
+      const mesh_region *rg;
       ga_tree *ptree;
       base_vector elem;
       tree_description(void) : ptree(0) {}
@@ -160,6 +160,13 @@
     };
 
   private:
+
+    std::map<const mesh *, std::list<mesh_region> > registred_mims;
+    mesh dummy_mesh;
+    mesh dummy_mim;
+    mesh_region dummy_region;
+    
+    mesh_region &register_region(const mesh &m, const mesh_region &region);
 
     typedef std::map<std::string, var_description> VAR_SET;
 

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4702&r1=4701&r2=4702&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Fri Jul  4 14:50:17 2014
@@ -2951,9 +2951,9 @@
   struct ga_instruction_copy_tensor_possibly_void : public ga_instruction {
     base_tensor &t, &tc1;
     virtual int exec(void) {
-      // static int count = 0;
       GA_DEBUG_INFO("Instruction: tensor copy possibly void");
       // cout << "copy tensor " << tc1 << endl;
+      // static int count = 0;
       // if (++count > 10) GMM_ASSERT1(false, "stop here");
       if (tc1.size())
         gmm::copy(tc1.as_vector(), t.as_vector());
@@ -3927,6 +3927,20 @@
   // functions, operators.
   //=========================================================================
 
+  mesh_region &ga_workspace::register_region(const mesh &m,
+                                             const mesh_region &region) {
+    if (&m == &dummy_mesh) return dummy_region;
+
+    std::list<mesh_region> &lmr = registred_mims[&m];
+    for (std::list<mesh_region>::iterator it = lmr.begin();
+         it != lmr.end(); ++it) {
+      if (it->compare(m, region, m)) return *it;
+    }
+    lmr.push_back(region);
+    return lmr.back();
+  }
+
+
   typedef std::pair<std::string, std::string> var_trans_pair;
   
   static bool ga_extract_variables(pga_tree_node pnode,
@@ -3987,6 +4001,7 @@
       //     ga_print_node(tree.root, cout); cout << endl;
       bool remain = true;
       size_type order = 0, ind_tree = 0;
+
       switch(tree.root->test_function_type) {
       case 0: order = 0; break;
       case 1: order = 1; break;
@@ -4010,7 +4025,7 @@
             trees[i].name_test2.compare(tree.root->name_test2) == 0 &&
             trees[i].interpolate_name_test2.compare
             (tree.root->interpolate_name_test2) == 0 &&
-            trees[i].rg.compare(*(trees[i].m), rg, m)) {
+            trees[i].rg == &rg) {
           ga_tree &ftree = *(trees[i].ptree);
             
           ftree.insert_node(ftree.root);
@@ -4028,7 +4043,7 @@
         ind_tree = trees.size(); remain = false;
         trees.push_back(tree_description());
         trees.back().mim = &mim; trees.back().m = &m;
-        trees.back().rg = rg;
+        trees.back().rg = &rg;
         trees.back().ptree = new ga_tree;
         trees.back().ptree->swap(tree);
         pga_tree_node root = trees.back().ptree->root;
@@ -4071,7 +4086,8 @@
 
   size_type ga_workspace::add_expression(const std::string expr,
                                          const mesh_im &mim,
-                                         const mesh_region &rg) {
+                                         const mesh_region &rg_) {
+    const mesh_region &rg = register_region(mim.linked_mesh(), rg_);
     // cout << "adding expression " << expr << endl;
     size_type max_order = 0;
     ga_tree tree;
@@ -4102,22 +4118,20 @@
   }
 
   void ga_workspace::add_function_expression(const std::string expr) {
-    static mesh_im dummy_mim;
-    static mesh dummy_mesh;
     ga_tree tree;
     ga_read_string(expr, tree);
     ga_semantic_analysis(expr, tree, *this, 1, false, true);
     if (tree.root) {
       // GMM_ASSERT1(tree.root->nb_test_functions() == 0,
       //            "Invalid function expression");
-      add_tree(tree, dummy_mesh, dummy_mim, 0, expr, false);
+      add_tree(tree, dummy_mesh, dummy_mim, dummy_region, expr, false);
     }
   }
 
   void ga_workspace::add_interpolation_expression(const std::string expr,
                                                   const mesh &m,
-                                                  mesh_region rg) {
-    static mesh_im dummy_mim;
+                                                  mesh_region rg_) {
+    const mesh_region &rg = register_region(m, rg_);
     ga_tree tree;
     ga_read_string(expr, tree);
     ga_semantic_analysis(expr, tree, *this, m.dim(), false, false);
@@ -7742,7 +7756,7 @@
       if (root) {
         GMM_ASSERT1(!scalar || (root->t.size() == 1),
                     "The result of the given expression is not a scalar"); 
-        ga_instruction_set::region_mim rm(td.mim, &(td.rg));
+        ga_instruction_set::region_mim rm(td.mim, td.rg);
         gis.whole_instructions[rm].m = td.m;
         ga_if_hierarchy if_hierarchy;
         ga_compile_node(root, workspace, gis,
@@ -7835,7 +7849,7 @@
         pga_tree_node root = gis.trees.back().root;
         if (root) {
           // Compile tree
-          ga_instruction_set::region_mim rm(td.mim, &(td.rg));
+          ga_instruction_set::region_mim rm(td.mim, td.rg);
           ga_instruction_set::region_mim_instructions &rmi
             = gis.whole_instructions[rm];
           rmi.m = td.m;
@@ -7872,7 +7886,7 @@
           // Compiling tree
           // cout << "compiling "; ga_print_node(root, cout); cout << endl;
 
-          ga_instruction_set::region_mim rm(td.mim, &(td.rg));
+          ga_instruction_set::region_mim rm(td.mim, td.rg);
           ga_instruction_set::region_mim_instructions &rmi
             = gis.whole_instructions[rm];
           rmi.m = td.m;
@@ -8123,6 +8137,7 @@
       
       const getfem::mesh_im &mim = *(it->first.first);
       const getfem::mesh &m = *(it->second.m);
+
       GMM_ASSERT1(&m == &(mim.linked_mesh()), "Incompatibility of meshes");
       size_type P = m.dim();
       ga_instruction_list &gil = it->second.instructions;