[Getfem-commits] r4993 - in /trunk/getfem: interface/tests/python/demo_Mindlin_Reissner_plate.py src/getfem_generic_assembly.cc

[Yves . Renard]

Author: renard
Date: Wed May 13 21:25:43 2015
New Revision: 4993

URL: http://svn.gna.org/viewcvs/getfem?rev=4993&view=rev
Log:
bug fix

Modified:
    trunk/getfem/interface/tests/python/demo_Mindlin_Reissner_plate.py
    trunk/getfem/src/getfem_generic_assembly.cc

Modified: 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=4993&r1=4992&r2=4993&view=diff
==============================================================================
--- trunk/getfem/interface/tests/python/demo_Mindlin_Reissner_plate.py  
(original)
+++ trunk/getfem/interface/tests/python/demo_Mindlin_Reissner_plate.py  Wed May 
13 21:25:43 2015
@@ -31,11 +31,11 @@
 import os
 
 ## Parameters
-Emodulus = 1           # Young Modulus
+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
+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
@@ -43,7 +43,7 @@
 K = 1                  # Degree of the finite element method
 dirichlet_version = 1  # 0 = simplification, 1 = with multipliers,
                        # 2 = penalization
-r = 1E8                # Penalization parameter.
+r = 1.E8               # Penalization parameter.
 NX = 80                # Number of element per direction
 
 
@@ -95,78 +95,3 @@
 
 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/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4993&r1=4992&r2=4993&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Wed May 13 21:25:43 2015
@@ -745,6 +745,32 @@
       break;
     case GA_NODE_CONSTANT: case GA_NODE_ZERO:
       if (pnode1->t.size() != pnode2->t.size()) return false;
+      
+      switch(version) {
+      case 0: case 1:
+        if (pnode1->test_function_type != pnode2->test_function_type)
+          return false;
+        if ((pnode1->test_function_type & 1) &&
+            pnode1->name_test1.compare(pnode2->name_test1) != 0)
+          return false;
+        if ((pnode1->test_function_type & 2) &&
+            pnode1->name_test2.compare(pnode2->name_test2) != 0)
+          return false;
+        break;
+      case 2:
+        if ((pnode1->test_function_type == 1 &&
+             pnode2->test_function_type == 1) ||
+            (pnode1->test_function_type == 2 &&
+             pnode2->test_function_type == 2))
+          return false;
+        if ((pnode1->test_function_type & 1) &&
+            pnode1->name_test1.compare(pnode2->name_test2) != 0)
+          return false;
+        if ((pnode1->test_function_type & 2) &&
+            pnode1->name_test2.compare(pnode2->name_test1) != 0)
+          return false;
+        break;
+      }
       if (pnode1->t.size() != 1 &&
           pnode1->t.sizes().size() != pnode2->t.sizes().size()) return false;
       for (size_type i = 0; i < pnode1->t.sizes().size(); ++i)
@@ -5081,7 +5107,12 @@
 
     switch (pnode->node_type) {
     case GA_NODE_CONSTANT: case GA_NODE_ZERO:
-      c += ga_hash_code(pnode->t); break;
+      c += ga_hash_code(pnode->t);
+      if (pnode->test_function_type & 1)
+        c += 34.731 * ga_hash_code(pnode->name_test1);
+      if (pnode->test_function_type & 2)
+        c += 34.731 * ga_hash_code(pnode->name_test2);
+      break;
 
     case GA_NODE_OP: c += scalar_type(pnode->op_type)*M_E*M_PI*M_PI; break;
     case GA_NODE_X: c += scalar_type(pnode->nbc1) + M_E*M_PI; break;
@@ -8215,10 +8246,11 @@
         // cout << " and "; ga_print_node(*it, cout); cout << endl;
         if (sub_tree_are_equal(pnode, *it, workspace, 1)) {
           // cout << "confirmed no transpose" << endl;
-          if (pnode->t.size() == 1)
+          if (pnode->t.size() == 1) {
             pgai = new ga_instruction_copy_scalar(pnode->t[0], (*it)->t[0]);
-          else
+          } else {
             pgai = new ga_instruction_copy_tensor(pnode->t, (*it)->t);
+          }
           rmi.instructions.push_back(pgai);
           return;
         }