[Getfem-commits] r4762 - in /trunk/getfem: doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst interface/src/gf_linsolve.cc

[logari81]

Author: logari81
Date: Wed Sep 10 11:51:18 2014
New Revision: 4762

URL: http://svn.gna.org/viewcvs/getfem?rev=4762&view=rev
Log:
add interface for the MUMPS solver in gf_linsolve

Modified:
    
trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst
    trunk/getfem/interface/src/gf_linsolve.cc

Modified: 
trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst?rev=4762&r1=4761&r2=4762&view=diff
==============================================================================
--- 
trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst 
    (original)
+++ 
trunk/getfem/doc/sphinx/source/userdoc/model_contact_friction_large_sliding.rst 
    Wed Sep 10 11:51:18 2014
@@ -413,7 +413,8 @@
   add_contact_boundary_to_large_sliding_contact_brick(model &md,
       size_type indbrick, const mesh_im &mim, size_type region,
       bool is_master, bool is_slave, const std::string &u,
-      const std::string &lambda = "", const std::string &w = "")
+      const std::string &lambda = "", const std::string &w = "",
+      bool frame_indifferent = false)
 
 where ``region`` should be a valid mesh region number representing a boundary, 
``is_master`` should be set to ``true`` if the contact detection is to be done 
on that contact boundary, ``is_slave`` should be set to ``true`` if the 
integration of contact terms is to be done on that boundary. Note that a 
contact boundary is allowed to be both master and slave, in particular to allow 
self-contact detection. ``u`` is the displacement variable. If ``is_slave`` is 
set to true, ``lambda`` should describe a multiplier variable with degrees of 
freedom on the contact boundary (typically added to the model with the 
``md.add_filtered_fem_variable(...) method). Pure master contact boundary do 
not need the definition of a multiplier. Additionally, ``w`` is for the 
evolutionnary case and represents the displacement at the previous time step.
 
@@ -423,5 +424,5 @@
       size_type indbrick, std::string expr, size_type N)
 
 where `expr` is an expression using the high-level
-generic assembly language (whith `X` is the current position) which should
-be a signed distance to the obstacle. `N` is the mesh dimension.
+generic assembly language (with `X` is the current position) which should
+be a signed distance to the obstacle. `N` is the mesh dimension.

Modified: trunk/getfem/interface/src/gf_linsolve.cc
URL: 
http://svn.gna.org/viewcvs/getfem/trunk/getfem/interface/src/gf_linsolve.cc?rev=4762&r1=4761&r2=4762&view=diff
==============================================================================
--- trunk/getfem/interface/src/gf_linsolve.cc   (original)
+++ trunk/getfem/interface/src/gf_linsolve.cc   Wed Sep 10 11:51:18 2014
@@ -25,6 +25,7 @@
 #include <gmm/gmm_iter_solvers.h>
 #include <getfemint_misc.h>
 #include <getfem/getfem_superlu.h>
+#include <gmm/gmm_MUMPS_interface.h>
 
 using namespace getfemint;
 
@@ -32,8 +33,8 @@
 
 template <typename T> static void
 iterative_gmm_solver(iterative_gmm_solver_type stype, gsparse &gsp,
-                   getfemint::mexargs_in& in,
-                    getfemint::mexargs_out& out, T) {
+                    getfemint::mexargs_in& in,
+                     getfemint::mexargs_out& out, T) {
   garray<T> b = in.pop().to_garray(int(gsp.nrows()), T());
   garray<T> x = out.pop().create_array_v(int(gsp.nrows()), T());
 
@@ -74,7 +75,7 @@
 }
 
 void iterative_gmm_solver(iterative_gmm_solver_type stype,
-                        getfemint::mexargs_in& in, getfemint::mexargs_out& 
out) {
+                         getfemint::mexargs_in& in, getfemint::mexargs_out& 
out) {
   dal::shared_ptr<gsparse> pgsp = in.pop().to_sparse();
   gsparse &gsp = *pgsp;
   if (!gsp.is_complex() && in.front().is_complex())
@@ -85,7 +86,7 @@
 
 template <typename T> static void
 superlu_solver(gsparse &gsp,
-              getfemint::mexargs_in& in, getfemint::mexargs_out& out, T) {
+               getfemint::mexargs_in& in, getfemint::mexargs_out& out, T) {
   garray<T> b = in.pop().to_garray(int(gsp.nrows()), T());
   garray<T> x = out.pop().create_array(b.getm(), b.getn(), T());
   double rcond;
@@ -95,6 +96,15 @@
     out.pop().from_scalar(rcond ? 1./rcond : 0.);
 }
 
+template <typename T> static void
+mumps_solver(gsparse &gsp,
+             getfemint::mexargs_in& in, getfemint::mexargs_out& out, T) {
+  garray<T> b = in.pop().to_garray(int(gsp.nrows()), T());
+  garray<T> x = out.pop().create_array(b.getm(), b.getn(), T());
+  gsp.to_csc();
+  gmm::MUMPS_solve(gsp.csc(T()),x,b);
+}
+
 /address@hidden
   Various linear system solvers.
 @*/
@@ -105,7 +115,7 @@
 struct sub_gf_linsolve : virtual public dal::static_stored_object {
   int arg_in_min, arg_in_max, arg_out_min, arg_out_max;
   virtual void run(getfemint::mexargs_in& in,
-                  getfemint::mexargs_out& out) = 0;
+                   getfemint::mexargs_out& out) = 0;
 };
 
 typedef boost::intrusive_ptr<sub_gf_linsolve> psub_command;
@@ -114,15 +124,15 @@
 template <typename T> static inline void dummy_func(T &) {}
 
 #define sub_command(name, arginmin, arginmax, argoutmin, argoutmax, code) { \
-    struct subc : public sub_gf_linsolve {                             \
-      virtual void run(getfemint::mexargs_in& in,                      \
-                      getfemint::mexargs_out& out)                     \
-      { dummy_func(in); dummy_func(out); code }                                
\
-    };                                                                 \
-    psub_command psubc = new subc;                                     \
-    psubc->arg_in_min = arginmin; psubc->arg_in_max = arginmax;                
\
-    psubc->arg_out_min = argoutmin; psubc->arg_out_max = argoutmax;    \
-    subc_tab[cmd_normalize(name)] = psubc;                             \
+    struct subc : public sub_gf_linsolve {                                  \
+      virtual void run(getfemint::mexargs_in& in,                           \
+                       getfemint::mexargs_out& out)                         \
+      { dummy_func(in); dummy_func(out); code }                             \
+    };                                                                      \
+    psub_command psubc = new subc;                                          \
+    psubc->arg_in_min = arginmin; psubc->arg_in_max = arginmax;             \
+    psubc->arg_out_min = argoutmin; psubc->arg_out_max = argoutmax;         \
+    subc_tab[cmd_normalize(name)] = psubc;                                  \
   }                           
 
 
@@ -168,11 +178,11 @@
     /address@hidden @CELL{U, cond} = ('lu', @tsp M, @vec b)
       Alias for ::LINSOLVE('superlu',...)@*/
     sub_command
-      ("lu", 2, 2, 0, 1,
+      ("lu", 2, 2, 0, 2,
        dal::shared_ptr<gsparse> pgsp = in.pop().to_sparse();
        gsparse &gsp = *pgsp;
        if (!gsp.is_complex() && in.front().is_complex())
-        THROW_BADARG("please use a real right hand side, or convert the sparse 
matrix to a complex one");
+         THROW_BADARG("please use a real right hand side, or convert the 
sparse matrix to a complex one");
        if (gsp.is_complex()) superlu_solver(gsp, in, out, complex_type());
        else                  superlu_solver(gsp, in, out, scalar_type());
        );
@@ -183,13 +193,25 @@
 
     The condition number estimate `cond` is returned with the solution 
address@hidden/
     sub_command
-      ("superlu", 2, 2, 0, 1,
+      ("superlu", 2, 2, 0, 2,
        dal::shared_ptr<gsparse> pgsp = in.pop().to_sparse();
        gsparse &gsp = *pgsp;
        if (!gsp.is_complex() && in.front().is_complex())
-        THROW_BADARG("please use a real right hand side, or convert the sparse 
matrix to a complex one");
+         THROW_BADARG("please use a real right hand side, or convert the 
sparse matrix to a complex one");
        if (gsp.is_complex()) superlu_solver(gsp, in, out, complex_type());
        else                  superlu_solver(gsp, in, out, scalar_type());
+       );
+
+    /address@hidden @CELL{U, cond} = ('mumps', @tsp M, @vec b)
+    Solve `M.U = b` using the MUMPS address@hidden/
+    sub_command
+      ("mumps", 2, 2, 0, 1,
+       dal::shared_ptr<gsparse> pgsp = in.pop().to_sparse();
+       gsparse &gsp = *pgsp;
+       if (!gsp.is_complex() && in.front().is_complex())
+         THROW_BADARG("please use a real right hand side, or convert the 
sparse matrix to a complex one");
+       if (gsp.is_complex()) mumps_solver(gsp, in, out, complex_type());
+       else                  mumps_solver(gsp, in, out, scalar_type());
        );
 
   }
@@ -205,8 +227,8 @@
   SUBC_TAB::iterator it = subc_tab.find(cmd);
   if (it != subc_tab.end()) {
     check_cmd(cmd, it->first.c_str(), m_in, m_out, it->second->arg_in_min,
-             it->second->arg_in_max, it->second->arg_out_min,
-             it->second->arg_out_max);
+              it->second->arg_in_max, it->second->arg_out_min,
+              it->second->arg_out_max);
     it->second->run(m_in, m_out);
   }
   else bad_cmd(init_cmd);