Home
Publications
Presentations
Support
Download
Examples
Developers
Class Docs
Namespaces
Classes
Files
Mailing Lists
Gallery
Wiki

condensed_eigen_system.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 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

#include "libmesh/libmesh_config.h"

// Currently, the EigenSystem should only be available
// if SLEPc support is enabled.
#if defined(LIBMESH_HAVE_SLEPC)

#include "libmesh/condensed_eigen_system.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/equation_systems.h"
#include "libmesh/dof_map.h"
#include "libmesh/parallel.h"

namespace libMesh
{

CondensedEigenSystem::CondensedEigenSystem (EquationSystems& es,
                                            const std::string& name,
                                            const unsigned int number)
  : Parent(es, name, number),
    condensed_matrix_A(SparseMatrix<Number>::build()),
    condensed_matrix_B(SparseMatrix<Number>::build()),
    condensed_dofs_initialized(false)
{
}

void CondensedEigenSystem::initialize_condensed_dofs(std::set<unsigned int>& global_dirichlet_dofs_set)
{
  // First, put all local dofs into non_dirichlet_dofs_set and
  std::set<unsigned int> local_non_condensed_dofs_set;
  for(unsigned int i=this->get_dof_map().first_dof(); i<this->get_dof_map().end_dof(); i++)
    local_non_condensed_dofs_set.insert(i);

  // Now erase the condensed dofs
  std::set<unsigned int>::iterator iter     = global_dirichlet_dofs_set.begin();
  std::set<unsigned int>::iterator iter_end = global_dirichlet_dofs_set.end();

  for ( ; iter != iter_end ; iter++)
  {
    unsigned int condensed_dof_index = *iter;
    if ( (this->get_dof_map().first_dof() <= condensed_dof_index) &&
         (condensed_dof_index < this->get_dof_map().end_dof()) )
    {
      local_non_condensed_dofs_set.erase(condensed_dof_index);
    }
  }

  // Finally, move local_non_condensed_dofs_set over to a vector for convenience in solve()
  iter     = local_non_condensed_dofs_set.begin();
  iter_end = local_non_condensed_dofs_set.end();

  this->local_non_condensed_dofs_vector.clear();

  for ( ; iter != iter_end; ++iter)
    {
      unsigned int non_condensed_dof_index = *iter;

      this->local_non_condensed_dofs_vector.push_back(non_condensed_dof_index);
    }

  condensed_dofs_initialized = true;
}

unsigned int CondensedEigenSystem::n_global_non_condensed_dofs() const
{
  if(!condensed_dofs_initialized)
  {
    return this->n_dofs();
  }
  else
  {
    unsigned int n_global_non_condensed_dofs = local_non_condensed_dofs_vector.size();
    CommWorld.sum(n_global_non_condensed_dofs);

    return n_global_non_condensed_dofs;
  }
}


void CondensedEigenSystem::solve()
{
  START_LOG("solve()", "CondensedEigenSystem");

  // If we haven't initialized any condensed dofs,
  // just use the default eigen_system
  if(!condensed_dofs_initialized)
  {
    STOP_LOG("solve()", "CondensedEigenSystem");
    Parent::solve();
    return;
  }

  // A reference to the EquationSystems
  EquationSystems& es = this->get_equation_systems();

  // check that necessary parameters have been set
  libmesh_assert (es.parameters.have_parameter<unsigned int>("eigenpairs"));
  libmesh_assert (es.parameters.have_parameter<unsigned int>("basis vectors"));

  if (this->assemble_before_solve)
    // Assemble the linear system
    this->assemble ();

  // If we reach here, then there should be some non-condensed dofs
  libmesh_assert(!local_non_condensed_dofs_vector.empty());

  // Now condense the matrices
  matrix_A->create_submatrix(*condensed_matrix_A,
                             local_non_condensed_dofs_vector,
                             local_non_condensed_dofs_vector);

  matrix_B->create_submatrix(*condensed_matrix_B,
                             local_non_condensed_dofs_vector,
                             local_non_condensed_dofs_vector);

  // Get the tolerance for the solver and the maximum
  // number of iterations. Here, we simply adopt the linear solver
  // specific parameters.
  const Real tol            =
    es.parameters.get<Real>("linear solver tolerance");

  const unsigned int maxits =
    es.parameters.get<unsigned int>("linear solver maximum iterations");

  const unsigned int nev    =
    es.parameters.get<unsigned int>("eigenpairs");

  const unsigned int ncv    =
    es.parameters.get<unsigned int>("basis vectors");

  std::pair<unsigned int, unsigned int> solve_data;

  // call the solver depending on the type of eigenproblem
  if ( generalized() )
    {
      //in case of a generalized eigenproblem
      solve_data = eigen_solver->solve_generalized
        (*condensed_matrix_A,*condensed_matrix_B, nev, ncv, tol, maxits);
    }

  else
    {
      libmesh_assert (!matrix_B);

      //in case of a standard eigenproblem
      solve_data = eigen_solver->solve_standard (*condensed_matrix_A, nev, ncv, tol, maxits);
    }

  set_n_converged(solve_data.first);
  set_n_iterations(solve_data.second);

  STOP_LOG("solve()", "CondensedEigenSystem");
}



std::pair<Real, Real> CondensedEigenSystem::get_eigenpair(unsigned int i)
{
  START_LOG("get_eigenpair()", "CondensedEigenSystem");

  // If we haven't initialized any condensed dofs,
  // just use the default eigen_system
  if(!condensed_dofs_initialized)
  {
    STOP_LOG("get_eigenpair()", "CondensedEigenSystem");
    return Parent::get_eigenpair(i);
  }

  // If we reach here, then there should be some non-condensed dofs
  libmesh_assert(!local_non_condensed_dofs_vector.empty());

  // This function assumes that condensed_solve has just been called.
  // If this is not the case, then we will trip an asset in get_eigenpair
  AutoPtr< NumericVector<Number> > temp = NumericVector<Number>::build();
  unsigned int n_local = local_non_condensed_dofs_vector.size();
  unsigned int n       = n_local;
  CommWorld.sum(n);

  temp->init (n, n_local, false, libMeshEnums::PARALLEL);

  std::pair<Real, Real> eval = eigen_solver->get_eigenpair (i, *temp);

  // Now map temp to solution. Loop over local entries of local_non_condensed_dofs_vector
  this->solution->zero();
  for(unsigned int i=0; i<local_non_condensed_dofs_vector.size(); i++)
  {
    unsigned int index = local_non_condensed_dofs_vector[i];
    solution->set(index,(*temp)(temp->first_local_index()+i));
  }

  solution->close();
  this->update();

  STOP_LOG("get_eigenpair()", "CondensedEigenSystem");

  return eval;
}




} // namespace libMesh


#endif // LIBMESH_HAVE_SLEPC

---

Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:45 UTC

Hosted By: