trilinos_epetra_vector.c

一个用来实现偏微分方程中网格的计算库

// $Id: epetra_vector.C 2606 2008-01-23 20:21:47Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// 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 "trilinos_epetra_vector.h"

#ifdef HAVE_TRILINOS

#include "parallel.h"
#include "utility.h"
#include "dense_vector.h"
#include "parallel.h"

// Trilinos Includes
#include <Epetra_Import.h>



template <typename T>
T EpetraVector<T>::sum () const
{
  libmesh_assert(this->closed());
  
  const unsigned int nl = _vec->MyLength();

  T sum=0.0;

  T * values = _vec->Values();
  
  for (unsigned int i=0; i<nl; i++)
    sum += values[i];

  Parallel::sum<T>(sum);
  
  return sum;
}

template <typename T>
Real EpetraVector<T>::l1_norm () const
{
  libmesh_assert(this->closed());

  Real value;

  _vec->Norm1(&value);
  
  return value;
}

template <typename T>
Real EpetraVector<T>::l2_norm () const
{
  libmesh_assert(this->closed());
  
  Real value;

  _vec->Norm2(&value);
  
  return value;
}

template <typename T>
Real EpetraVector<T>::linfty_norm () const
{
  libmesh_assert(this->closed());
  
  Real value;

  _vec->NormInf(&value);
  
  return value;
}

template <typename T>
NumericVector<T>&
EpetraVector<T>::operator += (const NumericVector<T>& v)
{
  libmesh_assert(this->closed());
  
  this->add(1., v);
  
  return *this;
}



template <typename T>
NumericVector<T>&
EpetraVector<T>::operator -= (const NumericVector<T>& v)
{
  libmesh_assert(this->closed());
  
  this->add(-1., v);
  
  return *this;
}



template <typename T>
void EpetraVector<T>::set (const unsigned int i_in, const T value_in)
{
  int i = static_cast<int> (i_in);
  T value = value_in;

  libmesh_assert(i_in<this->size());

  _vec->ReplaceGlobalValues(1, &i, &value);

  this->_is_closed = false;
}



template <typename T>
void EpetraVector<T>::add (const unsigned int i_in, const T value_in)
{
  int i = static_cast<int> (i_in);
  T value = value_in;

  libmesh_assert(i_in<this->size());
  
  _vec->SumIntoGlobalValues(1, &i, &value);

  this->_is_closed = false;
}



template <typename T>
void EpetraVector<T>::add_vector (const std::vector<T>& v,
				  const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (v.size() == dof_indices.size());

  _vec->SumIntoGlobalValues (v.size(),
			     (int*) &dof_indices[0],
			     const_cast<T*>(&v[0]));

}



template <typename T>
void EpetraVector<T>::add_vector (const NumericVector<T>& V,
				 const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (V.size() == dof_indices.size());

  for (unsigned int i=0; i<V.size(); i++)
    this->add (dof_indices[i], V(i));
}



// TODO: fill this in after creating an EpetraMatrix
template <typename T>
void EpetraVector<T>::add_vector (const NumericVector<T>& /* V_in */,
				  const SparseMatrix<T>& /* A_in */)
{
  LIBMESH_THROW(libMesh::NotImplemented());

//   const EpetraVector<T>* V = dynamic_cast<const EpetraVector<T>*>(&V_in);
//   const EpetraMatrix<T>* A = dynamic_cast<const EpetraMatrix<T>*>(&A_in);

//   libmesh_assert (V != NULL);
//   libmesh_assert (A != NULL);
  
//   int ierr=0;

//   A->close();

//   // The const_cast<> is not elegant, but it is required since Epetra
//   // is not const-correct.  
//   ierr = MatMultAdd(const_cast<EpetraMatrix<T>*>(A)->mat(), V->_vec, _vec, _vec);
//          CHKERRABORT(libMesh::COMM_WORLD,ierr); 
}



template <typename T>
void EpetraVector<T>::add_vector (const DenseVector<T>& /* V_in */,
				  const std::vector<unsigned int>& /* dof_indices */)
{
  LIBMESH_THROW(libMesh::NotImplemented());

//   libmesh_assert (V_in.size() == dof_indices.size());

//   const EpetraVector<T>* V = dynamic_cast<const EpetraVector<T>*>(&V_in);

//   this->add_vector(V->_vec->get_values[0], dof_indices);
}



template <typename T>
void EpetraVector<T>::add (const T v_in)
{
  const unsigned int nl = _vec->MyLength();

  T * values = _vec->Values();
  
  for (unsigned int i=0; i<nl; i++)
    values[i]+=v_in;

  this->_is_closed = false;
}


template <typename T>
void EpetraVector<T>::add (const NumericVector<T>& v)
{
  this->add (1., v);
}


template <typename T>
void EpetraVector<T>::add (const T a_in, const NumericVector<T>& v_in)
{
  const EpetraVector<T>* v = dynamic_cast<const EpetraVector<T>*>(&v_in);

  libmesh_assert (v != NULL);
  libmesh_assert(this->size() == v->size());

  _vec->Update(a_in,*v->_vec, 1.);
}



template <typename T>
void EpetraVector<T>::insert (const std::vector<T>& v,
			      const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (v.size() == dof_indices.size());

  _vec->ReplaceGlobalValues (v.size(),
			     (int*) &dof_indices[0],
			     const_cast<T*>(&v[0]));
}



template <typename T>
void EpetraVector<T>::insert (const NumericVector<T>& V,
			      const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (V.size() == dof_indices.size());

  // TODO: If V is an EpetraVector this can be optimized
  for (unsigned int i=0; i<V.size(); i++)
    this->set (dof_indices[i], V(i));
}



template <typename T>
void EpetraVector<T>::insert (const DenseVector<T>& v,
			      const std::vector<unsigned int>& dof_indices)
{
  libmesh_assert (v.size() == dof_indices.size());
  
  std::vector<T> &vals = const_cast<DenseVector<T>&>(v).get_values();
  
  _vec->ReplaceGlobalValues (v.size(),
			     (int*) &dof_indices[0],
			     &vals[0]);
}



template <typename T>
void EpetraVector<T>::scale (const T factor_in)
{
  _vec->Scale(factor_in);
}



template <typename T>
T EpetraVector<T>::dot (const NumericVector<T>& V_in) const
{
  const EpetraVector<T>* V = dynamic_cast<const EpetraVector<T>*>(&V_in);

  libmesh_assert(V);
  
  T result=0.0;

  _vec->Dot(*V->_vec, &result);

  return result;
}


template <typename T>
NumericVector<T>& 
EpetraVector<T>::operator = (const T s_in)
{
  _vec->PutScalar(s_in);
  
  return *this;
}



template <typename T>
NumericVector<T>&
EpetraVector<T>::operator = (const NumericVector<T>& v_in)
{
  const EpetraVector<T>* v = dynamic_cast<const EpetraVector<T>*>(&v_in);

  libmesh_assert (v != NULL);
  
  *this = *v;
  
  return *this;
}



template <typename T>
EpetraVector<T>&
EpetraVector<T>::operator = (const EpetraVector<T>& v)
{
  (*_vec) = *v._vec;
  
  return *this;
}



template <typename T>
NumericVector<T>&
EpetraVector<T>::operator = (const std::vector<T>& v)
{
  T * values = _vec->Values();
  
  /**
   * Case 1:  The vector is the same size of
   * The global vector.  Only add the local components.
   */
  if(this->size() == v.size())
  {
    const unsigned int nl=this->local_size();
    const unsigned int fli=this->first_local_index();