$extrastylesheet
petsc_auto_fieldsplit.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library.
00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
00003 
00004 // This library is free software; you can redistribute it and/or
00005 // modify it under the terms of the GNU Lesser General Public
00006 // License as published by the Free Software Foundation; either
00007 // version 2.1 of the License, or (at your option) any later version.
00008 
00009 // This library is distributed in the hope that it will be useful,
00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012 // Lesser General Public License for more details.
00013 
00014 // You should have received a copy of the GNU Lesser General Public
00015 // License along with this library; if not, write to the Free Software
00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00017 
00018 #include "libmesh/petsc_auto_fieldsplit.h"
00019 
00020 #ifdef LIBMESH_HAVE_PETSC
00021 
00022 EXTERN_C_FOR_PETSC_BEGIN
00023 #  include <petscksp.h>
00024 EXTERN_C_FOR_PETSC_END
00025 
00026 // Local includes
00027 #include "libmesh/dof_map.h"
00028 #include "libmesh/system.h"
00029 
00030 #if !PETSC_VERSION_LESS_THAN(3,2,0)
00031 
00032 // C++ includes
00033 
00034 namespace {
00035 using namespace libMesh;
00036 
00037 void indices_to_fieldsplit (const Parallel::Communicator& comm,
00038                             const std::vector<dof_id_type>& indices,
00039                             PC my_pc,
00040                             const std::string& field_name)
00041 {
00042   const PetscInt *idx = PETSC_NULL;
00043   if (!indices.empty())
00044     idx = reinterpret_cast<const PetscInt*>(&indices[0]);
00045 
00046   IS is;
00047   int ierr = ISCreateLibMesh(comm.get(), indices.size(),
00048                              idx, PETSC_COPY_VALUES, &is);
00049   CHKERRABORT(comm.get(), ierr);
00050 
00051   ierr = PCFieldSplitSetIS(my_pc, field_name.c_str(), is);
00052   CHKERRABORT(comm.get(), ierr);
00053 }
00054 
00055 }
00056 
00057 namespace libMesh
00058 {
00059 
00060 void petsc_auto_fieldsplit (PC my_pc,
00061                             const System &sys)
00062 {
00063   std::string sys_prefix = "--solver_group_";
00064 
00065   if (libMesh::on_command_line("--solver_system_names"))
00066     {
00067       sys_prefix = sys_prefix + sys.name() + "_";
00068     }
00069 
00070   std::map<std::string, std::vector<dof_id_type> > group_indices;
00071 
00072   if (libMesh::on_command_line("--solver_variable_names"))
00073     {
00074       for (unsigned int v = 0; v != sys.n_vars(); ++v)
00075         {
00076           const std::string& var_name = sys.variable_name(v);
00077 
00078           std::vector<dof_id_type> var_idx;
00079           sys.get_dof_map().local_variable_indices
00080             (var_idx, sys.get_mesh(), v);
00081 
00082           std::string group_command = sys_prefix + var_name;
00083 
00084           const std::string empty_string;
00085 
00086           std::string group_name = libMesh::command_line_value
00087             (group_command, empty_string);
00088 
00089           if (group_name != empty_string)
00090             {
00091               std::vector<dof_id_type> &indices =
00092                 group_indices[group_name];
00093               const bool prior_indices = !indices.empty();
00094               indices.insert(indices.end(), var_idx.begin(),
00095                              var_idx.end());
00096               if (prior_indices)
00097                 std::sort(indices.begin(), indices.end());
00098             }
00099           else
00100             {
00101               indices_to_fieldsplit (sys.comm(), var_idx, my_pc, var_name);
00102             }
00103         }
00104     }
00105 
00106   for (std::map<std::string, std::vector<dof_id_type> >::const_iterator
00107          i = group_indices.begin(); i != group_indices.end(); ++i)
00108     {
00109       indices_to_fieldsplit(sys.comm(), i->second, my_pc, i->first);
00110     }
00111 }
00112 
00113 } // namespace libMesh
00114 
00115 
00116 #else  // #PETSC_VERSION < 3.2.0
00117 
00118 namespace libMesh
00119 {
00120 void petsc_auto_fieldsplit (PC /* my_pc */,
00121                             const System & /* sys */)
00122 {
00123   if (libMesh::on_command_line("--solver_variable_names"))
00124     {
00125       libmesh_do_once(
00126                       libMesh::out << "WARNING: libMesh does not support setting field splits" <<
00127                       std::endl << "with PETSc "
00128                       << LIBMESH_DETECTED_PETSC_VERSION_MAJOR << '.'
00129                       << LIBMESH_DETECTED_PETSC_VERSION_MINOR << '.'
00130                       << LIBMESH_DETECTED_PETSC_VERSION_SUBMINOR << std::endl;);
00131     }
00132 }
00133 }
00134 
00135 #endif // #PETSC_VERSION > 3.2.0
00136 #endif // #ifdef LIBMESH_HAVE_PETSC