[Getfem-commits] r4901 - in /trunk/getfem/src: getfem_generic_assembly.cc getfem_linearized_plates.cc

[Yves . Renard]

Author: renard
Date: Wed Mar 25 11:13:06 2015
New Revision: 4901

URL: http://svn.gna.org/viewcvs/getfem?rev=4901&view=rev
Log:
working elementary projection but suboptimal

Modified:
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/src/getfem_linearized_plates.cc

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=4901&r1=4900&r2=4901&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Wed Mar 25 11:13:06 2015
@@ -319,8 +319,10 @@
     size_type pos;                // Position of the first character in string
     std::string name;             // variable/constant/function/operator name
     std::string interpolate_name; // For Interpolate : name of transformation
-    std::string interpolate_name_der; // For Interpolate frivative:
+    std::string interpolate_name_der; // For Interpolate derivative:
                                       // name of transformation
+    std::string elementary_name;  // For Elementary_transformation :
+                                  // name of transformation
     size_type der1, der2;         // For functions and nonlinear operators,
                                   // optional derivative or second derivative.
     GA_TOKEN_TYPE op_type;
@@ -754,6 +756,8 @@
     case GA_NODE_ELEMENTARY_GRAD_TEST: case GA_NODE_ELEMENTARY_HESS_TEST: 
       if (pnode1->interpolate_name.compare(pnode2->interpolate_name))
         return false;
+      if (pnode1->elementary_name.compare(pnode2->elementary_name))
+        return false;
       // The test continues with what follows
     case GA_NODE_TEST: case GA_NODE_GRAD_TEST: case GA_NODE_HESS_TEST:
       {
@@ -789,6 +793,8 @@
     case GA_NODE_INTERPOLATE_HESS: case GA_NODE_ELEMENTARY_VAL:
     case GA_NODE_ELEMENTARY_GRAD: case GA_NODE_ELEMENTARY_HESS: 
       if (pnode1->interpolate_name.compare(pnode2->interpolate_name))
+        return false;
+      if (pnode1->elementary_name.compare(pnode2->elementary_name))
         return false;
       if (pnode1->name.compare(pnode2->name)) return false;
       break;
@@ -995,7 +1001,7 @@
           << pnode->interpolate_name << ")"; break;
     case GA_NODE_ELEMENTARY:
       str << "Elementary_transformation(" << pnode->name << ","
-          << pnode->interpolate_name << ")"; break;
+          << pnode->elementary_name << ")"; break;
     case GA_NODE_INTERPOLATE_FILTER:
       str << "Interpolate_filter(" << pnode->interpolate_name << ",";
       ga_print_node(pnode->children[0], str);
@@ -1015,7 +1021,7 @@
       break;
     case GA_NODE_ELEMENTARY_VAL:
       str << "Elementary_transformation(" << pnode->name << ","
-          << pnode->interpolate_name << ")";
+          << pnode->elementary_name << ")";
       break;
     case GA_NODE_GRAD: str << "Grad_" << pnode->name; break;
     case GA_NODE_INTERPOLATE_GRAD:
@@ -1024,7 +1030,7 @@
       break;
     case GA_NODE_ELEMENTARY_GRAD:
       str << "Elementary_transformation(Grad_" << pnode->name << ","
-          << pnode->interpolate_name << ")";
+          << pnode->elementary_name << ")";
       break;
     case GA_NODE_HESS: str << "Hess_" << pnode->name; break;
     case GA_NODE_INTERPOLATE_HESS:
@@ -1033,7 +1039,7 @@
       break;
     case GA_NODE_ELEMENTARY_HESS:
       str << "Elementary_transformation(Hess_" << pnode->name << ","
-          << pnode->interpolate_name << ")";
+          << pnode->elementary_name << ")";
       break;
     case GA_NODE_TEST:
       if (pnode->test_function_type == 1) str << "Test_" << pnode->name;
@@ -1052,7 +1058,7 @@
     case GA_NODE_ELEMENTARY_TEST:
       str << "Elementary_transformation(";
       if (pnode->test_function_type == 1) str << "Test_"; else str << "Test2_";
-      str << pnode->name << "," << pnode->interpolate_name << ")";
+      str << pnode->name << "," << pnode->elementary_name << ")";
       break;
     case GA_NODE_GRAD_TEST:
       if (pnode->test_function_type == 1) str << "Grad_Test_" << pnode->name;
@@ -1066,7 +1072,7 @@
     case GA_NODE_ELEMENTARY_GRAD_TEST:
       str << "Elementary_transformation(Grad_";
       if (pnode->test_function_type == 1) str << "Test_"; else str << "Test2_";
-      str << pnode->name << "," << pnode->interpolate_name << ")";
+      str << pnode->name << "," << pnode->elementary_name << ")";
       break;
     case GA_NODE_HESS_TEST:
       if (pnode->test_function_type == 1) str << "Hess_Test_" << pnode->name;
@@ -1080,7 +1086,7 @@
     case GA_NODE_ELEMENTARY_HESS_TEST:
       str << "Elementary_transformation(Hess_";
       if (pnode->test_function_type == 1) str << "Test_"; else str << "Test2_";
-      str << pnode->name << "," << pnode->interpolate_name << ")";
+      str << pnode->name << "," << pnode->elementary_name << ")";
       break;
     case GA_NODE_SPEC_FUNC: str << pnode->name; break;
     case GA_NODE_OPERATOR:
@@ -1288,7 +1294,7 @@
               ga_throw_error(expr, pos,
                              "Second argument of Elementary_transformation 
should "
                              "be a transformation name.");
-            tree.current_node->interpolate_name
+            tree.current_node->elementary_name
               = std::string(&(expr[token_pos]), token_length);
             t_type = ga_get_token(expr, pos, token_pos, token_length);
             if (t_type != GA_RPAR)
@@ -3148,7 +3154,6 @@
       GA_DEBUG_ASSERT(t_.size() == Z.size() * Qmult * Qmult,
                       "Wrong size for base vector");
 
-
       t.adjust_sizes(t_.sizes()); 
 
       if (Qmult == 1) {
@@ -3172,8 +3177,11 @@
 
       gmm::resize(M, ndof*Qmult, ndof*Qmult);
       elemtrans->give_transformation(mf, ctx.convex_num(), M);
-      t.mat_transp_reduction(t_, M, 0);
-
+      // cout << "M = " << M << endl;
+      t_.mat_reduction(t, M, 0);
+      // cout << "t_ = " << t_ << endl;
+      // cout << "t = " << t << endl;
+      // gmm::copy(t.as_vector(), t_.as_vector());
       return 0;
     }
     ga_instruction_elementary_transformation_base
@@ -3225,7 +3233,7 @@
 
       gmm::resize(M, ndof*Qmult, ndof*Qmult);
       elemtrans->give_transformation(mf, ctx.convex_num(), M);
-      t.mat_transp_reduction(t_, M, 0); 
+      t_.mat_reduction(t, M, 0); 
 
       return 0;
     }
@@ -3279,7 +3287,7 @@
 
       gmm::resize(M, ndof*Qmult, ndof*Qmult);
       elemtrans->give_transformation(mf, ctx.convex_num(), M);
-      t.mat_transp_reduction(t_, M, 0);
+      t_.mat_reduction(t, M, 0);
 
       return 0;
     }
@@ -4412,21 +4420,16 @@
         pnode->node_type == GA_NODE_ELEMENTARY_HESS) {
       bool group = workspace.variable_group_exists(pnode->name);
       bool iscte = (!group) && workspace.is_constant(pnode->name);
-      const std::string &iname = (pnode->node_type == GA_NODE_INTERPOLATE_VAL 
||
-                                  pnode->node_type == GA_NODE_INTERPOLATE_GRAD 
||
-                                  pnode->node_type == GA_NODE_INTERPOLATE_HESS)
-        ? pnode->interpolate_name : std::string("");
-
       if (!iscte) found_var = true;
       if (!ignore_data || !iscte) {
         if (group && expand_groups) {
           const std::vector<std::string> &t
             = workspace.variable_group(pnode->name);
           for (size_type i = 0; i < t.size(); ++i)
-            vars.insert(var_trans_pair(t[i], iname));
+            vars.insert(var_trans_pair(t[i], pnode->interpolate_name));
 
         } else
-          vars.insert(var_trans_pair(pnode->name, iname));
+          vars.insert(var_trans_pair(pnode->name, pnode->interpolate_name));
       }
     }
     if (pnode->node_type == GA_NODE_INTERPOLATE_VAL ||
@@ -4940,7 +4943,8 @@
     case GA_NODE_ELEMENTARY_HESS: case GA_NODE_ELEMENTARY_TEST:
     case GA_NODE_ELEMENTARY_GRAD_TEST: case GA_NODE_ELEMENTARY_HESS_TEST:
       c += 1.33*(1.22+ga_hash_code(pnode->name))
-        + 1.66*ga_hash_code(pnode->interpolate_name);
+        + 1.66*ga_hash_code(pnode->interpolate_name)
+        + 2.63*ga_hash_code(pnode->elementary_name);
       break;
     case GA_NODE_INTERPOLATE_NORMAL: case GA_NODE_INTERPOLATE_X:
       c += M_PI*1.33*ga_hash_code(pnode->interpolate_name);
@@ -5382,7 +5386,7 @@
         }
 
         if (!(workspace.elementary_transformation_exists
-              (pnode->interpolate_name)))
+              (pnode->elementary_name)))
           ga_throw_error(expr, pnode->pos,
                          "Unknown elementary transformation");
       }
@@ -7974,21 +7978,21 @@
             pgai = new ga_instruction_elementary_transformation_val
               (pnode->t, rmi.base[mf],
                rmi.local_dofs[pnode->name], workspace.qdim(pnode->name),
-               workspace.elementary_transformation(pnode->interpolate_name),
+               workspace.elementary_transformation(pnode->elementary_name),
                *mf, gis.ctx);
             break;
           case GA_NODE_ELEMENTARY_GRAD:
             pgai = new ga_instruction_elementary_transformation_grad
               (pnode->t, rmi.grad[mf],
                rmi.local_dofs[pnode->name], workspace.qdim(pnode->name),
-               workspace.elementary_transformation(pnode->interpolate_name),
+               workspace.elementary_transformation(pnode->elementary_name),
                *mf, gis.ctx);
             break;
           case GA_NODE_ELEMENTARY_HESS:
             pgai = new ga_instruction_elementary_transformation_hess
               (pnode->t, rmi.hess[mf],
                rmi.local_dofs[pnode->name], workspace.qdim(pnode->name),
-               workspace.elementary_transformation(pnode->interpolate_name),
+               workspace.elementary_transformation(pnode->elementary_name),
                *mf, gis.ctx);
             break;
           default: break;
@@ -8102,19 +8106,19 @@
           case GA_NODE_ELEMENTARY_TEST:
             pgai = new ga_instruction_elementary_transformation_base
               (pnode->t, rmi.base[mf], mf->get_qdim(),
-               workspace.elementary_transformation(pnode->interpolate_name),
+               workspace.elementary_transformation(pnode->elementary_name),
                *mf, gis.ctx);
             break;
           case GA_NODE_ELEMENTARY_GRAD_TEST:
             pgai = new ga_instruction_elementary_transformation_grad_base
               (pnode->t, rmi.grad[mf], mf->get_qdim(),
-               workspace.elementary_transformation(pnode->interpolate_name),
+               workspace.elementary_transformation(pnode->elementary_name),
                *mf, gis.ctx);
             break;
           case GA_NODE_ELEMENTARY_HESS_TEST:
             pgai = new ga_instruction_elementary_transformation_hess_base
               (pnode->t, rmi.hess[mf], mf->get_qdim(),
-               workspace.elementary_transformation(pnode->interpolate_name),
+               workspace.elementary_transformation(pnode->elementary_name),
                *mf, gis.ctx);
             break;
           default: break;

Modified: trunk/getfem/src/getfem_linearized_plates.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_linearized_plates.cc?rev=4901&r1=4900&r2=4901&view=diff
==============================================================================
--- trunk/getfem/src/getfem_linearized_plates.cc        (original)
+++ trunk/getfem/src/getfem_linearized_plates.cc        Wed Mar 25 11:13:06 2015
@@ -86,13 +86,21 @@
         
       for (size_type i = 0; i < pim->nb_points_on_convex(); ++i) {
 
-        scalar_type coeff = pim->coeff(i); // Mult by ctx.J(), but not useful 
here
+        scalar_type coeff = pim->coeff(i); // Mult by ctx.J() not useful here
         ctx1.set_xref(pim->point(i));
         ctx2.set_xref(pim->point(i));    
         pf1->real_base_value(ctx1, t1);
         vectorize_base_tensor(t1, tv1, ndof1, pf1->target_dim(), N);
         pf2->real_base_value(ctx2, t2);
         vectorize_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
+
+        if (N == 2) { // H curl ... ! and for triangles ?
+          gmm::clear(tv2);
+          tv2(0,0) = scalar_type(1);
+          tv2(1,1) = scalar_type(1);
+          tv2(2,0) = pim->point(i)[1];
+          tv2(3,1) = pim->point(i)[0];
+        }
        
         gmm::mult(tv1, gmm::transposed(tv1), aux0);
         gmm::add(gmm::scaled(aux0, coeff), M1);
@@ -108,7 +116,8 @@
       gmm::lu_inverse(M2);
       gmm::mult(M1, B, aux1);
       gmm::mult(aux1, M2, aux2);
-      GMM_ASSERT1(gmm::mat_nrows(M) == ndof1, "Element not convenient for 
projection");
+      GMM_ASSERT1(gmm::mat_nrows(M) == ndof1,
+                  "Element not convenient for projection");
       gmm::mult(aux2, gmm::transposed(B), M);
       gmm::clean(M, 1E-15);
       // cout << "M = " << M << endl;