Dear Roy/Ben,
Please find attached the first version of "my part". I haven't tested
it because I had no idea how to do that. Well, at least it compiles.
Perhaps it is best that you guys implement your part now and then run
one of the examples on it, and in the case that it doesn't work as
expected, we will have to somehow cooperate in finding the bugs. If
you have a better idea, let me know.
Best Regards,
Tim
--
Dr. Tim Kroeger
[email protected] Phone +49-421-218-7710
[email protected] Fax +49-421-218-4236
Fraunhofer MEVIS, Institute for Medical Image Computing
Universitaetsallee 29, 28359 Bremen, Germany
Index: include/numerics/petsc_vector.h
===================================================================
--- include/numerics/petsc_vector.h (Revision 3206)
+++ include/numerics/petsc_vector.h (Arbeitskopie)
@@ -1,4 +1,4 @@
-// $Id$
+// $id: petsc_vector.h 3158 2008-11-17 17:44:57Z jwpeterson $
// The libMesh Finite Element Library.
// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
@@ -78,6 +78,15 @@
const unsigned int n_local);
/**
+ * Constructor. Set local dimension to \p n_local, the global
+ * dimension to \p n, but additionally reserve memory for the
+ * indices specified by the \p ghost argument.
+ */
+ PetscVector (const unsigned int N,
+ const unsigned int n_local,
+ const std::vector<unsigned int>& ghost);
+
+ /**
* Constructor. Creates a PetscVector assuming you already have a
* valid PETSc Vec object. In this case, v is NOT destroyed by the
* PetscVector constructor when this object goes out of scope.
@@ -136,6 +145,15 @@
void init (const unsigned int N,
const bool fast=false);
+ /**
+ * Create a vector that holds tha local indices plus those specified
+ * in the \p ghost argument.
+ */
+ virtual void init (const unsigned int /*N*/,
+ const unsigned int /*n_local*/,
+ const std::vector<unsigned int>& /*ghost*/,
+ const bool /*fast*/ = false);
+
// /**
// * Change the dimension to that of the
// * vector \p V. The same applies as for
@@ -233,6 +251,14 @@
* actually stored on this processor
*/
unsigned int last_local_index() const;
+
+ /**
+ * Maps the global index \p i to the corresponding global index. If
+ * the index is not a ghost cell, this is done by subtraction the
+ * number of the first local index. If it is a ghost cell, it has
+ * to be looked up in the map.
+ */
+ unsigned int map_global_to_local_index(const unsigned int i) const;
/**
* Access components, returns \p U(i).
@@ -441,7 +467,19 @@
*/
Vec _vec;
+
/**
+ * Type for map that maps global to local ghost cells.
+ */
+ typedef std::map<unsigned int,unsigned int> GlobalToLocalMap;
+
+ /**
+ * Map that maps global to local ghost cells (will be empty if not
+ * in ghost cell mode).
+ */
+ GlobalToLocalMap _global_to_local_map;
+
+ /**
* This boolean value should only be set to false
* for the constructor which takes a PETSc Vec object.
*/
@@ -456,7 +494,8 @@
template <typename T>
inline
PetscVector<T>::PetscVector ()
- : _destroy_vec_on_exit(true)
+ : _destroy_vec_on_exit(true),
+ _global_to_local_map()
{}
@@ -464,7 +503,8 @@
template <typename T>
inline
PetscVector<T>::PetscVector (const unsigned int n)
- : _destroy_vec_on_exit(true)
+ : _destroy_vec_on_exit(true),
+ _global_to_local_map()
{
this->init(n, n, false);
}
@@ -475,22 +515,60 @@
inline
PetscVector<T>::PetscVector (const unsigned int n,
const unsigned int n_local)
- : _destroy_vec_on_exit(true)
+ : _destroy_vec_on_exit(true),
+ _global_to_local_map()
{
this->init(n, n_local, false);
}
+template <typename T>
+inline
+PetscVector<T>::PetscVector (const unsigned int n,
+ const unsigned int n_local,
+ const std::vector<unsigned int>& ghost)
+ : _destroy_vec_on_exit(true),
+ _global_to_local_map()
+{
+ this->init(n, n_local, ghost, false);
+}
+
+
+
template <typename T>
inline
PetscVector<T>::PetscVector (Vec v)
- : _destroy_vec_on_exit(false)
+ : _destroy_vec_on_exit(false),
+ _global_to_local_map()
{
this->_vec = v;
this->_is_initialized = true;
+
+ /* We need to ask PETSc about the (local to global) ghost value
+ mapping and create the inverse mapping out of it. */
+ int ierr=0, petsc_size=0, petsc_local_size=0;
+ ierr = VecGetSize(_vec, &petsc_size);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGetLocalSize(_vec, &petsc_local_size);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ /* \p petsc_local_size is the number of non-ghost values. If it
+ equals the global size, then we are a serial vector, and there
+ are no ghost values. */
+ if(petsc_size!=petsc_local_size)
+ {
+ const unsigned int local_size = static_cast<unsigned
int>(petsc_local_size);
+ ISLocalToGlobalMapping mapping = _vec->mapping;
+ const unsigned int ghost_begin = static_cast<unsigned
int>(petsc_local_size);
+ const unsigned int ghost_end = static_cast<unsigned int>(mapping->n);
+ for(unsigned int i=ghost_begin; i<ghost_end; i++)
+ {
+ _global_to_local_map[mapping->indices[i]] = i-local_size;
+ }
+ }
}
@@ -564,6 +642,46 @@
template <typename T>
inline
+void PetscVector<T>::init (const unsigned int n,
+ const unsigned int n_local,
+ const std::vector<unsigned int>& ghost,
+ const bool fast)
+{
+ int ierr=0;
+ int petsc_n=static_cast<int>(n);
+ int petsc_n_local=static_cast<int>(n_local);
+ int petsc_n_ghost=static_cast<int>(ghost.size());
+ int* petsc_ghost = const_cast<int*>(reinterpret_cast<const int*>(&ghost[0]));
+
+ // Clear initialized vectors
+ if (this->initialized())
+ this->clear();
+
+ /* Make the global-to-local ghost cell map. */
+ for(unsigned int i=0; i<ghost.size(); i++)
+ {
+ _global_to_local_map[ghost[i]] = i;
+ }
+
+ /* Create vector. */
+ ierr = VecCreateGhost (libMesh::COMM_WORLD, petsc_n_local, petsc_n,
+ petsc_n_ghost, petsc_ghost,
+ &_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ ierr = VecSetFromOptions (_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+
+ this->_is_initialized = true;
+
+ if (fast == false)
+ this->zero ();
+}
+
+
+
+template <typename T>
+inline
void PetscVector<T>::close ()
{
libmesh_assert (this->initialized());
@@ -571,10 +689,18 @@
int ierr=0;
ierr = VecAssemblyBegin(_vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
ierr = VecAssemblyEnd(_vec);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ if(_global_to_local_map.size()!=0)
+ {
+ ierr = VecGhostUpdateBegin(_vec,ADD_VALUES,SCATTER_REVERSE);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostUpdateEnd(_vec,ADD_VALUES,SCATTER_REVERSE);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+
this->_is_closed = true;
}
@@ -593,6 +719,8 @@
}
this->_is_closed = this->_is_initialized = false;
+
+ _global_to_local_map.clear();
}
@@ -706,24 +834,61 @@
template <typename T>
inline
+unsigned int PetscVector<T>::map_global_to_local_index (const unsigned int i)
const
+{
+ libmesh_assert (this->initialized());
+
+ const unsigned int first = this->first_local_index();
+ const unsigned int last = this->last_local_index();
+
+ if((i>=first) && (i<last))
+ {
+ return i-first;
+ }
+
+ GlobalToLocalMap::const_iterator it = _global_to_local_map.find(i);
+ libmesh_assert (it!=_global_to_local_map.end());
+ return it->second+last-first;
+}
+
+
+
+template <typename T>
+inline
T PetscVector<T>::operator() (const unsigned int i) const
{
+ const unsigned int local_index = this->map_global_to_local_index(i);
libmesh_assert (this->initialized());
- libmesh_assert ( ((i >= this->first_local_index()) &&
- (i < this->last_local_index())) );
int ierr=0;
- PetscScalar *values, value=0.;
-
+ PetscScalar value=0.;
- ierr = VecGetArray(_vec, &values);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ if(_global_to_local_map.empty())
+ {
+ PetscScalar *values;
+ ierr = VecGetArray(_vec, &values);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ value = values[local_index];
+ ierr = VecRestoreArray (_vec, &values);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
+ else
+ {
+ /* Vectors that include ghost values require a special
+ handling. */
+ Vec loc_vec;
+ PetscScalar *values;
+ ierr = VecGhostGetLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGetArray(loc_vec, &values);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ value = values[local_index];
+ ierr = VecRestoreArray (loc_vec, &values);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ ierr = VecGhostRestoreLocalForm (_vec,&loc_vec);
+ CHKERRABORT(libMesh::COMM_WORLD,ierr);
+ }
- value = values[i - this->first_local_index()];
-
- ierr = VecRestoreArray (_vec, &values);
- CHKERRABORT(libMesh::COMM_WORLD,ierr);
-
return static_cast<T>(value);
}
@@ -771,6 +936,7 @@
{
std::swap(_vec, v._vec);
std::swap(_destroy_vec_on_exit, v._destroy_vec_on_exit);
+ std::swap(_global_to_local_map, v._global_to_local_map);
}
Index: include/numerics/numeric_vector.h
===================================================================
--- include/numerics/numeric_vector.h (Revision 3206)
+++ include/numerics/numeric_vector.h (Arbeitskopie)
@@ -71,6 +71,15 @@
*/
NumericVector (const unsigned n,
const unsigned int n_local);
+
+ /**
+ * Constructor. Set local dimension to \p n_local, the global
+ * dimension to \p n, but additionally reserve memory for the
+ * indices specified by the \p ghost argument.
+ */
+ NumericVector (const unsigned int N,
+ const unsigned int n_local,
+ const std::vector<unsigned int>& ghost);
public:
@@ -144,6 +153,15 @@
virtual void init (const unsigned int,
const bool = false) {}
+ /**
+ * Create a vector that holds tha local indices plus those specified
+ * in the \p ghost argument.
+ */
+ virtual void init (const unsigned int /*N*/,
+ const unsigned int /*n_local*/,
+ const std::vector<unsigned int>& /*ghost*/,
+ const bool /*fast*/ = false) {}
+
// /**
// * Change the dimension to that of the
// * vector \p V. The same applies as for
@@ -252,7 +270,8 @@
/**
* @returns the local size of the vector
- * (index_stop-index_start)
+ * (index_stop-index_start).
+ * In ghost cell mode, this does *not* include the ghost cells.
*/
virtual unsigned int local_size() const = 0;
@@ -556,6 +575,19 @@
template <typename T>
inline
+NumericVector<T>::NumericVector (const unsigned int n,
+ const unsigned int n_local,
+ const std::vector<unsigned int>& ghost) :
+ _is_closed(false),
+ _is_initialized(false)
+{
+ init(n, n_local, ghost, false);
+}
+
+
+
+template <typename T>
+inline
NumericVector<T>::~NumericVector ()
{
clear ();
Index: src/numerics/petsc_vector.C
===================================================================
--- src/numerics/petsc_vector.C (Revision 3206)
+++ src/numerics/petsc_vector.C (Arbeitskopie)
@@ -244,7 +244,7 @@
int ig = fli + i;
- PetscScalar value = (values[ig] + v);
+ PetscScalar value = (values[i] + v);
ierr = VecRestoreArray (_vec, &values);
CHKERRABORT(libMesh::COMM_WORLD,ierr);
@@ -441,6 +441,7 @@
if (v.initialized())
{
this->init (v.size(), v.local_size());
+ this->_global_to_local_map = v._global_to_local_map;
this->_is_closed = v._is_closed;
this->_is_initialized = v._is_initialized;
------------------------------------------------------------------------------
This SF.net email is sponsored by:
SourcForge Community
SourceForge wants to tell your story.
http://p.sf.net/sfu/sf-spreadtheword
_______________________________________________
Libmesh-devel mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-devel
