[Getfem-commits] r5302 - in /trunk/getfem: interface/src/ src/ src/getfem/

[Yves . Renard]

Author: renard
Date: Thu Apr 21 10:03:14 2016
New Revision: 5302

URL: http://svn.gna.org/viewcvs/getfem?rev=5302&view=rev
Log:
New method for a local projection of a scalar quantity onto a discontinuous 
mesh_fem (Von Mises stress for plasticity for instance). Could be improved

Modified:
    trunk/getfem/interface/src/gf_model_get.cc
    trunk/getfem/src/getfem/getfem_generic_assembly.h
    trunk/getfem/src/getfem_generic_assembly.cc
    trunk/getfem/src/getfem_models.cc

Modified: trunk/getfem/interface/src/gf_model_get.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_model_get.cc?rev=5302&r1=5301&r2=5302&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_model_get.cc  (original)
+++ trunk/getfem/interface/src/gf_model_get.cc  Thu Apr 21 10:03:14 2016
@@ -273,6 +273,26 @@
          getfem::ga_interpolation_mti(*md, expr, mti,
                                       result, rg, extrapolation, rg_source);
        }
+       out.pop().from_dcvector(result);
+       );
+
+    /address@hidden V = ('local_projection', @tmim mim, @str expr, @tmf mf[, 
@int region])
+      Make an elementwise L2 projection of an expression with respect
+      to the mesh_fem `mf`. This mesh_fem has to be
+      a discontinuous one.
+      The expression has to be valid according to the high-level generic
+      assembly language possibly including references to the variables
+      and data of the model. @*/
+    sub_command
+      ("local_projection", 3, 4, 0, 1,
+       getfem::mesh_im *mim = to_meshim_object(in.pop());
+       std::string expr = in.pop().to_string();
+       const getfem::mesh_fem *mf = to_meshfem_object(in.pop());
+       GMM_ASSERT1(!(mf->is_reduced()), "Sorry, cannot apply to reduced fems");
+       size_type region = in.remaining() ? in.pop().to_integer():size_type(-1);
+
+       getfem::base_vector result;
+       getfem::ga_local_projection(*md, *mim, expr, *mf, result, region);
        out.pop().from_dcvector(result);
        );
     

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=5302&r1=5301&r2=5302&view=diff
==============================================================================
--- trunk/getfem/src/getfem/getfem_generic_assembly.h   (original)
+++ trunk/getfem/src/getfem/getfem_generic_assembly.h   Thu Apr 21 10:03:14 2016
@@ -1,7 +1,7 @@
 /* -*- c++ -*- (enables emacs c++ mode) */
 /*===========================================================================
 
- Copyright (C) 2013-2015 Yves Renard
+ Copyright (C) 2013-2016 Yves Renard
 
  This file is a part of GetFEM++
 
@@ -429,11 +429,9 @@
     virtual ~ga_interpolation_context() {}
   };
 
-
   //=========================================================================
   // Interpolation functions
   //=========================================================================
-
 
   void ga_interpolation(ga_workspace &workspace,
                         ga_interpolation_context &gic);
@@ -448,7 +446,8 @@
   void ga_interpolation_mti
   (const getfem::model &md, const std::string &expr, mesh_trans_inv &mti,
    base_vector &result, const mesh_region &rg=mesh_region::all_convexes(),
-   int extrapolation = 0, const mesh_region 
&rg_source=mesh_region::all_convexes(),
+   int extrapolation = 0,
+   const mesh_region &rg_source=mesh_region::all_convexes(),
    size_type nbdof_ = size_type(-1));
 
   void ga_interpolation_im_data
@@ -458,6 +457,21 @@
   void ga_interpolation_im_data
   (ga_workspace &workspace, const im_data &imd, base_vector &result,
    const mesh_region &rg=mesh_region::all_convexes());
+
+  //=========================================================================
+  // Local projection functions
+  //=========================================================================
+
+  /** Make an elementwise L2 projection of an expression with respect
+      to the mesh_fem `mf`. This mesh_fem has to be a discontinuous one.
+      The expression has to be valid according to the high-level generic
+      assembly language possibly including references to the variables
+      and data of the model. 
+  */
+  void ga_local_projection(const getfem::model &md, const mesh_im &mim,
+                          const std::string &expr, const mesh_fem &mf,
+                          base_vector &result,
+                          const mesh_region &rg=mesh_region::all_convexes());
 
   //=========================================================================
   // Interpolate transformations

Modified: trunk/getfem/src/getfem_generic_assembly.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_generic_assembly.cc?rev=5302&r1=5301&r2=5302&view=diff
==============================================================================
--- trunk/getfem/src/getfem_generic_assembly.cc (original)
+++ trunk/getfem/src/getfem_generic_assembly.cc Thu Apr 21 10:03:14 2016
@@ -7358,9 +7358,10 @@
         break;
 
       case GA_DIV:
-        if (child1->tensor_order() > 0)
+        if (child1->tensor_proper_size() > 1)
           ga_throw_error(expr, pnode->pos,
-                         "Only the division by a scalar is allowed.");
+                         "Only the division by a scalar is allowed. "
+                        "Got a size of " << child1->tensor_proper_size());
         if (child1->test_function_type)
           ga_throw_error(expr, pnode->pos,
                          "Division by test functions is not allowed.");
@@ -7836,7 +7837,9 @@
         size_type s2 = (nbargs == 2) ? child2->t.size() : s1;
         if (s1 != s2 && (s1 != 1 || s2 != 1))
           ga_throw_error(expr, pnode->pos,
-                         "Invalid argument size for a scalar function.");
+                         "Invalid argument size for a scalar function. "
+                        "Size of first argument: " << s1 <<
+                        ". Size of second argument: " << s2 << ".");
 
         if (nbargs == 1) {
           pnode->t = child1->t;
@@ -11592,6 +11595,50 @@
     ga_interpolation(workspace, gic);
   }
 
+  //=========================================================================
+  // Local projection functions
+  //=========================================================================
+
+  void ga_local_projection(const getfem::model &md, const mesh_im &mim,
+                          const std::string &expr, const mesh_fem &mf,
+                          base_vector &result, const mesh_region &region) {
+
+    // The case where the expression is a vector one and mf a scalar fem is
+    // not taken into account for the moment.
+
+    // Could be improved by not make the assembly of the global mass matrix
+    // working locally. This means a specific assembly.
+    model_real_sparse_matrix M(mf.nb_dof(), mf.nb_dof());
+    asm_mass_matrix(M, mim, mf, region);
+       
+    ga_workspace workspace(md);
+    size_type nbdof = md.nb_dof();
+    gmm::sub_interval I(nbdof, mf.nb_dof());
+    workspace.add_fem_variable("c__dummy_var_95_", mf, I, base_vector(nbdof));
+    workspace.add_expression("("+expr+").Test_c__dummy_var_95_",
+                            mim, region, 2);
+    base_vector residual(nbdof+mf.nb_dof());
+    workspace.set_assembled_vector(residual);
+    workspace.assembly(1);
+    getfem::base_vector F(mf.nb_dof());
+    gmm::resize(result, mf.nb_dof());
+    gmm::copy(gmm::sub_vector(residual, I), F);
+
+    mesh_region rg(region);
+    mf.linked_mesh().intersect_with_mpi_region(rg);
+
+    getfem::base_matrix loc_M;
+    getfem::base_vector loc_U;
+    for (mr_visitor v(rg); !v.finished(); v.next()) {
+      size_type nd =  mf.nb_basic_dof_of_element(v.cv());
+      gmm::resize(loc_M, nd, nd); gmm::resize(loc_U, nd);
+      gmm::sub_index J(mf.ind_basic_dof_of_element(v.cv()));
+      gmm::copy(gmm::sub_matrix(M, J, J), loc_M);
+      gmm::lu_solve(loc_M, loc_U, gmm::sub_vector(F, J));
+      gmm::copy(loc_U, gmm::sub_vector(result, J));
+    }
+    MPI_SUM_VECTOR(result);
+  }
 
   //=========================================================================
   // Interpolate transformation with an expression

Modified: trunk/getfem/src/getfem_models.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/src/getfem_models.cc?rev=5302&r1=5301&r2=5302&view=diff
==============================================================================
--- trunk/getfem/src/getfem_models.cc   (original)
+++ trunk/getfem/src/getfem_models.cc   Thu Apr 21 10:03:14 2016
@@ -2741,7 +2741,6 @@
 
     // Generic expressions
     if (generic_expressions.size()) {
-
       model_real_plain_vector residual;
       if (version & BUILD_RHS) gmm::resize(residual, gmm::vect_size(rrhs));
 
@@ -2757,7 +2756,11 @@
         {
           exception.run([&]
           {
-            GMM_TRACE2("Global generic assembly");
+           if (version & BUILD_RHS)
+             GMM_TRACE2("Global generic assembly RHS");
+           if (version & BUILD_MATRIX)
+             GMM_TRACE2("Global generic assembly tangent term");
+
             ga_workspace workspace(*this);
 
             for (const auto &ge : generic_expressions)