Here is the current state now. I have no more questions (other than
the one which I asked in the previous mail). If Jed says that
everything looks well now, I will do the thing for the missing
overloaded PetscLinearSolver::solve() methods.
--
Dr. Tim Kroeger
CeVis -- Center of Complex Systems and Visualization
University of Bremen [email protected]
Universitaetsallee 29 [email protected]
D-28359 Bremen Phone +49-421-218-7710
Germany Fax +49-421-218-4236
Index: include/enums/enum_subset_solve_mode.h
===================================================================
--- include/enums/enum_subset_solve_mode.h (Revision 0)
+++ include/enums/enum_subset_solve_mode.h (Revision 0)
@@ -0,0 +1,45 @@
+// $Id: enum_norm_type.h 3430 2009-07-17 15:34:59Z roystgnr $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+#ifndef __enum_subset_solve_mode_h__
+#define __enum_subset_solve_mode_h__
+
+// ------------------------------------------------------------
+// enum SubsetSolveMode definition
+namespace libMeshEnums {
+
+ /**
+ * \enum libMeshEnums::SubsetSolveMode defines an \p enum for the
+ * question what happens to the dofs outside the given subset when a
+ * system is solved on a subset.
+ */
+ enum SubsetSolveMode {
+ SUBSET_ZERO = 0, //!< Set dofs outside the subset to zero.
+ SUBSET_COPY_RHS, //!< Set dofs outside the subset to the value of the
corresponding dofs of the right hand side.
+ SUBSET_DONT_TOUCH //!< Leaves dofs outside the subset unchanged. This is
fastest, but also most confusing because it abandons the property that the
solution vector is (theoretically) independent of the initial guess.
+ };
+
+}
+
+using namespace libMeshEnums;
+
+#endif // #ifndef __subset_solve_mode_h__
+
Index: include/systems/system_subset_by_subdomain.h
===================================================================
--- include/systems/system_subset_by_subdomain.h (Revision 0)
+++ include/systems/system_subset_by_subdomain.h (Revision 0)
@@ -0,0 +1,174 @@
+// $Id: system.h 3943 2010-08-26 13:28:04Z knezed01 $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+#ifndef __system_subset_by_subdomain_h__
+#define __system_subset_by_subdomain_h__
+
+// C++ includes
+#include <set>
+
+// Local Includes
+#include "system_subset.h"
+#include "id_types.h"
+
+namespace libMesh
+{
+
+// Forward Declarations
+
+/**
+ * This class represents a subset of the dofs of a \p System, selected
+ * by the \p subdomain_id and possible the variable numbers. The dofs
+ * in the subset will be sorted.
+ *
+ * @author Tim Kroeger, 2010.
+ */
+
+// ------------------------------------------------------------
+// SystemSubset class definition
+class SystemSubsetBySubdomain : public SystemSubset
+{
+public:
+
+ /**
+ * Subclass for user-specified selection of subdomain ids to be
+ * included in a \p SystemSubset.
+ */
+ class SubdomainSelection : public ReferenceCountedObject<SubdomainSelection>
+ {
+ public:
+
+ /**
+ * Constructor.
+ */
+ SubdomainSelection (void);
+
+ /**
+ * Destructor.
+ */
+ virtual ~SubdomainSelection (void);
+
+ /**
+ * Method that decides whether a given subdomain id is included in
+ * the subset or nor.
+ */
+ virtual bool operator()(const subdomain_id_type& subdomain_id)const=0;
+
+ private:
+ /**
+ * This isn't a copyable object, so let's make sure nobody tries.
+ *
+ * We won't even bother implementing this; we'll just make sure
+ * that the compiler doesn't implement a default.
+ */
+ SubdomainSelection(const SubdomainSelection&);
+
+ /**
+ * This isn't a copyable object, so let's make sure nobody tries.
+ *
+ * We won't even bother implementing this; we'll just make sure
+ * that the compiler doesn't implement a default.
+ */
+ SubdomainSelection& operator=(const SubdomainSelection&);
+
+ }; // subclass \p SubdomainSelection
+
+ /**
+ * Selection of subdomain ids by a list.
+ */
+ class SubdomainSelectionByList : public SubdomainSelection
+ {
+ public:
+ /**
+ * Constructor. Does not take a copy of the \p list, so make sure
+ * that the \p list does not go out of scope before \p *this does.
+ */
+ SubdomainSelectionByList (const std::set<subdomain_id_type>& list);
+
+ /**
+ * Method that decides whether a given subdomain id is included in
+ * the subset or nor.
+ */
+ virtual bool operator()(const subdomain_id_type& subdomain_id)const;
+
+ protected:
+ /**
+ * The actual list.
+ */
+ const std::set<subdomain_id_type>& _list;
+ };
+
+ /**
+ * Constructor. The subset will consist of those dofs which are
+ * associated to at least one mesh element that has a subdomain id
+ * contained in the \p subdomain_selection. If \p var_nums is not a
+ * \p NULL pointer, dofs that are associated to a variable number
+ * that is not contained in \p var_nums will not contain to the
+ * subset, no matter what elements they belong to.
+ */
+ SystemSubsetBySubdomain (const System& system,
+ const SubdomainSelection& subdomain_selection,
+ const std::set<unsigned int>* const var_nums = NULL);
+
+ /**
+ * Constructor. The subset will consist of those dofs which are
+ * associated to at least one mesh element that has a subdomain id
+ * contained in the set \p subdomain_ids. If \p var_nums is not a
+ * \p NULL pointer, dofs that are associated to a variable number
+ * that is not contained in \p var_nums will not contain to the
+ * subset, no matter what elements they belong to.
+ */
+ SystemSubsetBySubdomain (const System& system,
+ const std::set<subdomain_id_type>& subdomain_ids,
+ const std::set<unsigned int>* const var_nums = NULL);
+
+ /**
+ * Destructor.
+ */
+ virtual ~SystemSubsetBySubdomain ();
+
+ /**
+ * Method that returns the actual set of dofs that the subset
+ * consists of. The result will contain local dofs on each
+ * processor only and will not contain duplictates.
+ */
+ virtual const std::vector<unsigned int>& dof_ids(void)const;
+
+protected:
+
+ /**
+ * Initializes the class. Will be called by the constructors.
+ */
+ void init (const System& system,
+ const SubdomainSelection& subdomain_selection,
+ const std::set<unsigned int>* const var_nums);
+
+ /**
+ * The actual set of the dof ids.
+ */
+ std::vector<unsigned int> _dof_ids;
+
+}; // class SystemSubset
+
+} // namespace libMesh
+
+#endif
+
Index: include/systems/system.h
===================================================================
--- include/systems/system.h (Revision 3957)
+++ include/systems/system.h (Arbeitskopie)
@@ -31,6 +31,7 @@
#include "elem_range.h"
#include "enum_norm_type.h"
#include "enum_xdr_mode.h"
+#include "enum_subset_solve_mode.h"
#include "fe_type.h"
#include "libmesh_common.h"
#include "tensor_value.h" // For point_hessian
@@ -56,6 +57,7 @@
template <typename T> class VectorValue;
typedef VectorValue<Number> NumberVectorValue;
typedef NumberVectorValue Gradient;
+class SystemSubset;
/**
* This is the base class for classes which contain
@@ -162,6 +164,14 @@
* This method is only implemented in some derived classes.
*/
virtual void assemble_residual_derivatives (const ParameterVector&
parameters);
+
+ /**
+ * After calling this method, any solve will be restricted to the
+ * given subdomain. To disable this mode, call this method with \p
+ * subset being a \p NULL pointer.
+ */
+ virtual void restrict_solve_to (const SystemSubset* subset,
+ const SubsetSolveMode
subset_solve_mode=SUBSET_ZERO);
/**
* Solves the system. Must be overloaded in derived systems.
Index: include/systems/linear_implicit_system.h
===================================================================
--- include/systems/linear_implicit_system.h (Revision 3957)
+++ include/systems/linear_implicit_system.h (Arbeitskopie)
@@ -100,6 +100,14 @@
virtual void assemble () { ImplicitSystem::assemble(); }
/**
+ * After calling this method, any solve will be limited to the given
+ * subset. To disable this mode, call this method with \p subset
+ * being a \p NULL pointer.
+ */
+ virtual void restrict_solve_to (const SystemSubset* subset,
+ const SubsetSolveMode
subset_solve_mode=SUBSET_ZERO);
+
+ /**
* Assembles & solves the linear system A*x=b.
*/
virtual void solve ();
@@ -187,6 +195,16 @@
*/
ShellMatrix<Number>* _shell_matrix;
+ /**
+ * The current subset on which to solve (or \p NULL if none).
+ */
+ const SystemSubset* _subset;
+
+ /**
+ * If restrict-solve-to-subset mode is active, this member decides
+ * what happens with the dofs outside the subset.
+ */
+ SubsetSolveMode _subset_solve_mode;
};
Index: include/systems/system_subset.h
===================================================================
--- include/systems/system_subset.h (Revision 0)
+++ include/systems/system_subset.h (Revision 0)
@@ -0,0 +1,103 @@
+// $Id: system.h 3943 2010-08-26 13:28:04Z knezed01 $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+#ifndef __system_subset_h__
+#define __system_subset_h__
+
+// C++ includes
+#include <vector>
+
+// Local Includes
+#include "reference_counted_object.h"
+
+namespace libMesh
+{
+
+// Forward Declarations
+class System;
+
+/**
+ * This is a base class for classes which represent subsets of the
+ * dofs of a \p System.
+ *
+ * @author Tim Kroeger, 2010.
+ */
+
+// ------------------------------------------------------------
+// SystemSubset class definition
+class SystemSubset : public ReferenceCountedObject<SystemSubset>
+{
+public:
+
+ /**
+ * Constructor.
+ */
+ SystemSubset (const System& system);
+
+public:
+
+ /**
+ * Destructor.
+ */
+ virtual ~SystemSubset (void);
+
+ /**
+ * Method that returns the actual set of dofs that the subset
+ * consists of. The result must contain local dofs on each
+ * processor only and must not contain duplictates.
+ */
+ virtual const std::vector<unsigned int>& dof_ids(void)const=0;
+
+ /**
+ * Returns the \p System to which we belong.
+ */
+ const System& get_system(void)const;
+
+protected:
+
+ /**
+ * A reference to the \p System we belong to.
+ */
+ const System& _system;
+
+private:
+ /**
+ * This isn't a copyable object, so let's make sure nobody tries.
+ *
+ * We won't even bother implementing this; we'll just make sure that
+ * the compiler doesn't implement a default.
+ */
+ SystemSubset(const SystemSubset&);
+
+ /**
+ * This isn't a copyable object, so let's make sure nobody tries.
+ *
+ * We won't even bother implementing this; we'll just make sure that
+ * the compiler doesn't implement a default.
+ */
+ SystemSubset& operator=(const SystemSubset&);
+
+}; // class SystemSubset
+
+} // namespace libMesh
+
+#endif
+
Index: include/solvers/petsc_linear_solver.h
===================================================================
--- include/solvers/petsc_linear_solver.h (Revision 3957)
+++ include/solvers/petsc_linear_solver.h (Arbeitskopie)
@@ -122,7 +122,18 @@
* Initialize data structures if not done so already plus much more
*/
void init (PetscMatrix<T>* matrix);
+
/**
+ * After calling this method, all successive solves will be
+ * restricted to the given set of dofs, which must contain local
+ * dofs on each processor only and not contain any duplicates. This
+ * mode can be disabled by calling this method with \p dofs being a
+ * \p NULL pointer.
+ */
+ virtual void restrict_solve_to (const std::vector<unsigned int>* const dofs,
+ const SubsetSolveMode
subset_solve_mode=SUBSET_ZERO);
+
+ /**
* Call the Petsc solver. It calls the method below, using the
* same matrix for the system and preconditioner matrices.
*/
@@ -253,13 +264,34 @@
* Krylov subspace context
*/
KSP _ksp;
+
+ /**
+ * PETSc index set containing the dofs on which to solve (\p NULL
+ * means solve on all dofs).
+ */
+ IS _restrict_solve_to_is;
+
+ /**
+ * Internal method that returns the local size of \p
+ * _restrict_solve_to_is.
+ */
+ size_t _restrict_solve_to_is_local_size(void)const;
+
+ /**
+ * If restrict-solve-to-subset mode is active, this member decides
+ * what happens with the dofs outside the subset.
+ */
+ SubsetSolveMode _subset_solve_mode;
+
};
/*----------------------- functions ----------------------------------*/
template <typename T>
inline
-PetscLinearSolver<T>::PetscLinearSolver ()
+ PetscLinearSolver<T>::PetscLinearSolver ():
+ _restrict_solve_to_is(NULL),
+ _subset_solve_mode(SUBSET_ZERO)
{
if (libMesh::n_processors() == 1)
this->_preconditioner_type = ILU_PRECOND;
@@ -277,6 +309,23 @@
}
+
+template <typename T>
+inline size_t
+PetscLinearSolver<T>::
+_restrict_solve_to_is_local_size(void)const
+{
+ libmesh_assert(_restrict_solve_to_is!=NULL);
+
+ PetscInt s;
+ int ierr = ISGetLocalSize(_restrict_solve_to_is,&s);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ return static_cast<size_t>(s);
+}
+
+
+
} // namespace libMesh
Index: include/solvers/linear_solver.h
===================================================================
--- include/solvers/linear_solver.h (Revision 3957)
+++ include/solvers/linear_solver.h (Arbeitskopie)
@@ -24,12 +24,14 @@
// C++ includes
+#include <vector>
// Local includes
#include "libmesh_common.h"
#include "enum_solver_package.h"
#include "enum_solver_type.h"
#include "enum_preconditioner_type.h"
+#include "enum_subset_solve_mode.h"
#include "reference_counted_object.h"
#include "libmesh.h"
@@ -117,6 +119,16 @@
void attach_preconditioner(Preconditioner<T> * preconditioner);
/**
+ * After calling this method, all successive solves will be
+ * restricted to the given set of dofs, which must contain local
+ * dofs on each processor only and not contain any duplicates. This
+ * mode can be disabled by calling this method with \p dofs being a
+ * \p NULL pointer.
+ */
+ virtual void restrict_solve_to (const std::vector<unsigned int>* const dofs,
+ const SubsetSolveMode
subset_solve_mode=SUBSET_ZERO);
+
+ /**
* This function calls the solver
* "_solver_type" preconditioned with the
* "_preconditioner_type" preconditioner. Note that this method
Index: src/systems/linear_implicit_system.C
===================================================================
--- src/systems/linear_implicit_system.C (Revision 3957)
+++ src/systems/linear_implicit_system.C (Arbeitskopie)
@@ -28,6 +28,7 @@
#include "numeric_vector.h" // for parameter sensitivity calcs
#include "parameter_vector.h"
#include "sparse_matrix.h" // for get_transpose
+#include "system_subset.h"
namespace libMesh
{
@@ -43,7 +44,9 @@
linear_solver (LinearSolver<Number>::build()),
_n_linear_iterations (0),
_final_linear_residual (1.e20),
- _shell_matrix(NULL)
+ _shell_matrix(NULL),
+ _subset(NULL),
+ _subset_solve_mode(SUBSET_ZERO)
{
}
@@ -61,6 +64,8 @@
{
// clear the linear solver
linear_solver->clear();
+
+ this->restrict_solve_to(NULL);
// clear the parent data
Parent::clear();
@@ -79,6 +84,19 @@
+void LinearImplicitSystem::restrict_solve_to (const SystemSubset* subset,
+ const SubsetSolveMode
subset_solve_mode)
+{
+ _subset = subset;
+ _subset_solve_mode = subset_solve_mode;
+ if(subset!=NULL)
+ {
+ libmesh_assert(&subset->get_system()==this);
+ }
+}
+
+
+
void LinearImplicitSystem::solve ()
{
if (this->assemble_before_solve)
@@ -101,6 +119,11 @@
const unsigned int maxits =
es.parameters.get<unsigned int>("linear solver maximum iterations");
+ if(_subset!=NULL)
+ {
+ linear_solver->restrict_solve_to(&_subset->dof_ids(),_subset_solve_mode);
+ }
+
// Solve the linear system. Several cases:
std::pair<unsigned int, Real> rval = std::make_pair(0,0.0);
if(_shell_matrix)
@@ -110,6 +133,11 @@
// 2.) No shell matrix, with or without user-supplied preconditioner
rval = linear_solver->solve (*matrix,
this->request_matrix("Preconditioner"), *solution, *rhs, tol, maxits);
+ if(_subset!=NULL)
+ {
+ linear_solver->restrict_solve_to(NULL);
+ }
+
// Store the number of linear iterations required to
// solve and the final residual.
_n_linear_iterations = rval.first;
Index: src/systems/system_subset.C
===================================================================
--- src/systems/system_subset.C (Revision 0)
+++ src/systems/system_subset.C (Revision 0)
@@ -0,0 +1,52 @@
+// $Id: system.C 3943 2010-08-26 13:28:04Z knezed01 $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+// C++ includes
+
+
+// Local includes
+#include "system_subset.h"
+
+namespace libMesh
+{
+
+// ------------------------------------------------------------
+// SystemSubset implementation
+ SystemSubset::SystemSubset (const System& system):
+ _system(system)
+ {
+ }
+
+
+ SystemSubset::~SystemSubset (void)
+ {
+ }
+
+
+ const System&
+ SystemSubset::get_system(void)const
+ {
+ return _system;
+ }
+
+
+} // namespace libMesh
+
Index: src/systems/system_subset_by_subdomain.C
===================================================================
--- src/systems/system_subset_by_subdomain.C (Revision 0)
+++ src/systems/system_subset_by_subdomain.C (Revision 0)
@@ -0,0 +1,148 @@
+// $Id: system.C 3943 2010-08-26 13:28:04Z knezed01 $
+
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+
+// C++ includes
+
+
+// Local includes
+#include "system_subset_by_subdomain.h"
+#include "system.h"
+#include "dof_map.h"
+
+namespace libMesh
+{
+// ------------------------------------------------------------
+// SubdomainSelectio implementation
+ SystemSubsetBySubdomain::SubdomainSelection::
+ SubdomainSelection (void)
+ {
+ }
+
+ SystemSubsetBySubdomain::SubdomainSelection::
+ ~SubdomainSelection (void)
+ {
+ }
+
+ SystemSubsetBySubdomain::SubdomainSelectionByList::
+ SubdomainSelectionByList (const std::set<subdomain_id_type>& list):
+ _list(list)
+ {
+ }
+
+ bool
+ SystemSubsetBySubdomain::SubdomainSelectionByList::
+ operator()(const subdomain_id_type& subdomain_id)const
+ {
+ return _list.find(subdomain_id)!=_list.end();
+ }
+
+// ------------------------------------------------------------
+// SystemSubsetBySubdomain implementation
+
+ SystemSubsetBySubdomain::
+ SystemSubsetBySubdomain (const System& system,
+ const SubdomainSelection& subdomain_selection,
+ const std::set<unsigned int>* const var_nums):
+ SystemSubset(system),
+ _dof_ids()
+ {
+ this->init(system,subdomain_selection,var_nums);
+ }
+
+ SystemSubsetBySubdomain::
+ SystemSubsetBySubdomain (const System& system,
+ const std::set<subdomain_id_type>& subdomain_ids,
+ const std::set<unsigned int>* const var_nums):
+ SystemSubset(system),
+ _dof_ids()
+ {
+ SubdomainSelectionByList selection(subdomain_ids);
+ this->init(system,selection,var_nums);
+ }
+
+ SystemSubsetBySubdomain::
+ ~SystemSubsetBySubdomain (void)
+ {
+ }
+
+ const std::vector<unsigned int>&
+ SystemSubsetBySubdomain::
+ dof_ids(void)const
+ {
+ return _dof_ids;
+ }
+
+ void
+ SystemSubsetBySubdomain::
+ init (const System& system,
+ const SubdomainSelection& subdomain_selection,
+ const std::set<unsigned int>* const var_nums)
+ {
+ std::set<unsigned int> effective_var_nums;
+ if(var_nums!=NULL)
+ {
+ effective_var_nums = *var_nums;
+ }
+ else
+ {
+ for(unsigned int i=0; i<system.n_vars(); i++)
+ {
+ effective_var_nums.insert(i);
+ }
+ }
+
+ const DofMap & dof_map = system.get_dof_map();
+ std::vector<unsigned int> dof_indices;
+
+ const MeshBase& mesh = system.get_mesh();
+ MeshBase::const_element_iterator el =
mesh.active_local_elements_begin();
+ const MeshBase::const_element_iterator end_el =
mesh.active_local_elements_end();
+ for ( ; el != end_el; ++el)
+ {
+ const Elem* elem = *el;
+ if(subdomain_selection(elem->subdomain_id()))
+ {
+ std::set<unsigned int>::const_iterator it =
effective_var_nums.begin();
+ const std::set<unsigned int>::const_iterator itEnd =
effective_var_nums.end();
+ for (; it!=itEnd; ++it)
+ {
+ dof_map.dof_indices (elem, dof_indices, *it);
+ for(size_t i=0; i<dof_indices.size(); i++)
+ {
+ if((dof_indices[i]>=dof_map.variable_first_local_dof(*it))
&&
+ (dof_indices[i]<dof_map.variable_last_local_dof(*it)))
+ {
+ _dof_ids.push_back(dof_indices[i]);
+ }
+ }
+ }
+ }
+ }
+
+ /* Sort and unique the vector (using the same mechanism as in \p
+ DofMap::prepare_send_list()). */
+ std::sort(_dof_ids.begin(), _dof_ids.end());
+ std::vector<unsigned int>::iterator new_end = std::unique
(_dof_ids.begin(), _dof_ids.end());
+ std::vector<unsigned int> (_dof_ids.begin(), new_end).swap (_dof_ids);
+ }
+
+} // namespace libMesh
+
Index: src/systems/system.C
===================================================================
--- src/systems/system.C (Revision 3957)
+++ src/systems/system.C (Arbeitskopie)
@@ -367,6 +367,17 @@
+void System::restrict_solve_to (const SystemSubset* subset,
+ const SubsetSolveMode /*subset_solve_mode*/)
+{
+ if(subset!=NULL)
+ {
+ libmesh_not_implemented();
+ }
+}
+
+
+
void System::assemble ()
{
// Log how long the user's matrix assembly code takes
Index: src/solvers/petsc_linear_solver.C
===================================================================
--- src/solvers/petsc_linear_solver.C (Revision 3957)
+++ src/solvers/petsc_linear_solver.C (Arbeitskopie)
@@ -98,6 +98,8 @@
{
if (this->initialized())
{
+ this->restrict_solve_to(NULL);
+
this->_is_initialized = false;
int ierr=0;
@@ -294,7 +296,7 @@
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set operators. The input matrix works as the preconditioning matrix
- ierr = KSPSetOperators(_ksp, matrix->mat(),
matrix->mat(),SAME_NONZERO_PATTERN);
+ ierr = KSPSetOperators(_ksp, matrix->mat(),
matrix->mat(),DIFFERENT_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set user-specified solver and preconditioner types
@@ -360,9 +362,42 @@
+template <typename T>
+void
+PetscLinearSolver<T>::restrict_solve_to (const std::vector<unsigned int>*
const dofs,
+ const SubsetSolveMode
subset_solve_mode)
+{
+ libmesh_assert(this->initialized());
+ int ierr=0;
+ if(_restrict_solve_to_is!=NULL)
+ {
+ ierr = ISDestroy(_restrict_solve_to_is);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ _restrict_solve_to_is = NULL;
+ }
+ _subset_solve_mode = subset_solve_mode;
+
+ if(dofs!=NULL)
+ {
+ PetscInt* petsc_dofs;
+ ierr = PetscMalloc(dofs->size()*sizeof(PetscInt),&petsc_dofs);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ for(size_t i=0; i<dofs->size(); i++)
+ {
+ petsc_dofs[i] = (*dofs)[i];
+ }
+
+ ierr =
ISCreateGeneralNC(libMesh::COMM_WORLD,dofs->size(),petsc_dofs,&_restrict_solve_to_is);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+}
+
+
+
template <typename T>
std::pair<unsigned int, Real>
PetscLinearSolver<T>::solve (SparseMatrix<T>& matrix_in,
@@ -406,10 +441,15 @@
// 2.1.x & earlier style
#if PETSC_VERSION_LESS_THAN(2,2,0)
+
+ if(_restrict_solve_to_is!=NULL)
+ {
+ libmesh_not_implemented();
+ }
// Set operators. The input matrix works as the preconditioning matrix
ierr = SLESSetOperators(_sles, matrix->mat(), precond->mat(),
- SAME_NONZERO_PATTERN);
+ DIFFERENT_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the tolerances for the iterative solver. Use the user-supplied
@@ -431,6 +471,11 @@
// 2.2.0
#elif PETSC_VERSION_LESS_THAN(2,2,1)
+ if(_restrict_solve_to_is!=NULL)
+ {
+ libmesh_not_implemented();
+ }
+
// Set operators. The input matrix works as the preconditioning matrix
//ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
// SAME_NONZERO_PATTERN);
@@ -475,19 +520,65 @@
// 2.2.1 & newer style
#else
- // Set operators. The input matrix works as the preconditioning matrix
- if(!this->same_preconditioner)
- {
- ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
- SAME_NONZERO_PATTERN);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
- }
+ Mat submat = NULL;
+ Mat subprecond = NULL;
+ Vec subrhs;
+ Vec subsolution;
+ VecScatter scatter;
+
+ // Set operators. Also restrict rhs and solution vector to
+ // subdomain if neccessary.
+ if(_restrict_solve_to_is!=NULL)
+ {
+ size_t is_local_size = this->_restrict_solve_to_is_local_size();
+
+ ierr = VecCreate(libMesh::COMM_WORLD,&subrhs);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecSetFromOptions(subrhs);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecCreate(libMesh::COMM_WORLD,&subsolution);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecSetFromOptions(subsolution);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs,NULL,
&scatter);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr =
VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr =
VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr =
VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr =
VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = MatGetSubMatrix(matrix->mat(),
+ _restrict_solve_to_is,_restrict_solve_to_is,
+ PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = MatGetSubMatrix(precond->mat(),
+ _restrict_solve_to_is,_restrict_solve_to_is,
+ PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = KSPSetOperators(_ksp, submat, subprecond,
+ this->same_preconditioner ? SAME_PRECONDITIONER :
DIFFERENT_NONZERO_PATTERN);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
else
- {
- ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
- SAME_PRECONDITIONER);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
- }
+ {
+ ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
+ this->same_preconditioner ? SAME_PRECONDITIONER :
DIFFERENT_NONZERO_PATTERN);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
// Set the tolerances for the iterative solver. Use the user-supplied
// tolerance for the relative residual & leave the others at default values.
@@ -496,8 +587,16 @@
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Solve the linear system
- ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ if(_restrict_solve_to_is!=NULL)
+ {
+ ierr = KSPSolve (_ksp, subrhs, subsolution);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+ else
+ {
+ ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
// Get the number of iterations required for convergence
ierr = KSPGetIterationNumber (_ksp, &its);
@@ -507,6 +606,43 @@
ierr = KSPGetResidualNorm (_ksp, &final_resid);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ if(_restrict_solve_to_is!=NULL)
+ {
+ switch(_subset_solve_mode)
+ {
+ case SUBSET_ZERO:
+ ierr = VecZeroEntries(solution->vec());
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ break;
+
+ case SUBSET_COPY_RHS:
+ ierr = VecCopy(rhs->vec(),solution->vec());
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ break;
+
+ case SUBSET_DONT_TOUCH:
+ /* Nothing to do here. */
+ break;
+
+ }
+ ierr =
VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr =
VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecScatterDestroy(scatter);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecDestroy(subsolution);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecDestroy(subrhs);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = MatDestroy(submat);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = MatDestroy(subprecond);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+
#endif
STOP_LOG("solve()", "PetscLinearSolver");
@@ -574,7 +710,7 @@
// Set operators. The input matrix works as the preconditioning matrix
ierr = KSPSetOperators(_ksp, mat, mat,
- SAME_NONZERO_PATTERN);
+ DIFFERENT_NONZERO_PATTERN);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
// Set the tolerances for the iterative solver. Use the user-supplied
Index: src/solvers/linear_solver.C
===================================================================
--- src/solvers/linear_solver.C (Revision 3957)
+++ src/solvers/linear_solver.C (Arbeitskopie)
@@ -115,7 +115,19 @@
_preconditioner = preconditioner;
}
+template <typename T>
+void
+LinearSolver<T>::restrict_solve_to(const std::vector<unsigned int>* const dofs,
+ const SubsetSolveMode /*subset_solve_mode*/)
+{
+ if(dofs!=NULL)
+ {
+ libmesh_not_implemented();
+ }
+}
+
+
//------------------------------------------------------------------
// Explicit instantiations
template class LinearSolver<Number>;
------------------------------------------------------------------------------
Start uncovering the many advantages of virtual appliances
and start using them to simplify application deployment and
accelerate your shift to cloud computing.
http://p.sf.net/sfu/novell-sfdev2dev
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
