Epetra_FEVector.cpp

Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright (2001) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ************************************************************************
00028 //@HEADER
00029 
00030 #include <Epetra_FEVector.h>
00031 
00032 #include <Epetra_LocalMap.h>
00033 #include <Epetra_Comm.h>
00034 #include <Epetra_Map.h>
00035 #include <Epetra_Import.h>
00036 #include <Epetra_Export.h>
00037 #include <Epetra_Util.h>
00038 #include <Epetra_IntSerialDenseVector.h>
00039 #include <Epetra_SerialDenseVector.h>
00040 
00041 //----------------------------------------------------------------------------
00042 Epetra_FEVector::Epetra_FEVector(const Epetra_BlockMap& Map,
00043                                  int numVectors,
00044          bool ignoreNonLocalEntries)
00045   : Epetra_MultiVector(Map, numVectors),
00046     myFirstID_(0),
00047     myNumIDs_(0),
00048     nonlocalIDs_(NULL),
00049     nonlocalElementSize_(NULL),
00050     numNonlocalIDs_(0),
00051     numNonlocalIDsAlloc_(0),
00052     nonlocalCoefs_(NULL),
00053     numNonlocalCoefs_(0),
00054     numNonlocalCoefsAlloc_(0),
00055     ignoreNonLocalEntries_(ignoreNonLocalEntries)
00056 {
00057   myFirstID_ = Map.MinMyGID();
00058   myNumIDs_ = Map.NumMyElements();
00059   nonlocalCoefs_ = new double*[numVectors];
00060   for(int i=0; i<numVectors; ++i) nonlocalCoefs_[i] = NULL;
00061 }
00062 
00063 //----------------------------------------------------------------------------
00064 Epetra_FEVector::Epetra_FEVector(const Epetra_FEVector& source)
00065   : Epetra_MultiVector(source),
00066     myFirstID_(0),
00067     myNumIDs_(0),
00068     nonlocalIDs_(NULL),
00069     nonlocalElementSize_(NULL),
00070     numNonlocalIDs_(0),
00071     numNonlocalIDsAlloc_(0),
00072     nonlocalCoefs_(NULL),
00073     numNonlocalCoefsAlloc_(0),
00074     ignoreNonLocalEntries_(source.ignoreNonLocalEntries_)
00075 {
00076   *this = source;
00077 }
00078 
00079 //----------------------------------------------------------------------------
00080 Epetra_FEVector::~Epetra_FEVector()
00081 {
00082   destroyNonlocalData();
00083 
00084   delete [] nonlocalCoefs_;
00085   nonlocalCoefs_ = NULL;
00086 }
00087 
00088 //----------------------------------------------------------------------------
00089 int Epetra_FEVector::SumIntoGlobalValues(int numIDs, const int* GIDs,
00090                        const double* values,
00091                                          int vectorIndex)
00092 {
00093   return( inputValues( numIDs, GIDs, values, true, vectorIndex) );
00094 }
00095 
00096 //----------------------------------------------------------------------------
00097 int Epetra_FEVector::SumIntoGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
00098                        const Epetra_SerialDenseVector& values,
00099                                          int vectorIndex)
00100 {
00101   if (GIDs.Length() != values.Length()) {
00102     return(-1);
00103   }
00104 
00105   return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), true,
00106                        vectorIndex ) );
00107 }
00108 
00109 //----------------------------------------------------------------------------
00110 int Epetra_FEVector::SumIntoGlobalValues(int numIDs, const int* GIDs,
00111            const int* numValuesPerID,
00112                        const double* values,
00113                                          int vectorIndex)
00114 {
00115   return( inputValues( numIDs, GIDs, numValuesPerID, values, true,
00116                        vectorIndex) );
00117 }
00118 
00119 //----------------------------------------------------------------------------
00120 int Epetra_FEVector::ReplaceGlobalValues(int numIDs, const int* GIDs,
00121                        const double* values,
00122                                          int vectorIndex)
00123 {
00124   return( inputValues( numIDs, GIDs, values, false,
00125                        vectorIndex) );
00126 }
00127 
00128 //----------------------------------------------------------------------------
00129 int Epetra_FEVector::ReplaceGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
00130                        const Epetra_SerialDenseVector& values,
00131                                          int vectorIndex)
00132 {
00133   if (GIDs.Length() != values.Length()) {
00134     return(-1);
00135   }
00136 
00137   return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), false,
00138                        vectorIndex) );
00139 }
00140 
00141 //----------------------------------------------------------------------------
00142 int Epetra_FEVector::inputValues(int numIDs,
00143                                  const int* GIDs,
00144                                  const double* values,
00145                                  bool suminto,
00146                                  int vectorIndex)
00147 {
00148  //Important note!! This method assumes that there is only 1 point
00149  //associated with each element (GID), and writes to offset 0 in that
00150  //GID's block.
00151 
00152   for(int i=0; i<numIDs; ++i) {
00153     if (Map().MyGID(GIDs[i])) {
00154       if (suminto) {
00155         SumIntoGlobalValue(GIDs[i], 0, vectorIndex, values[i]);
00156       }
00157       else {
00158         ReplaceGlobalValue(GIDs[i], 0, vectorIndex, values[i]);
00159       }
00160     }
00161     else {
00162       if (!ignoreNonLocalEntries_) {
00163   EPETRA_CHK_ERR( inputNonlocalValue(GIDs[i], values[i], suminto,
00164                                            vectorIndex) );
00165       }
00166     }
00167   }
00168 
00169   return(0);
00170 }
00171 
00172 //----------------------------------------------------------------------------
00173 int Epetra_FEVector::ReplaceGlobalValues(int numIDs, const int* GIDs,
00174            const int* numValuesPerID,
00175                        const double* values,
00176                                          int vectorIndex)
00177 {
00178   return( inputValues( numIDs, GIDs, numValuesPerID, values, false,
00179                        vectorIndex) );
00180 }
00181 
00182 //----------------------------------------------------------------------------
00183 int Epetra_FEVector::inputValues(int numIDs,
00184                                  const int* GIDs,
00185          const int* numValuesPerID,
00186                                  const double* values,
00187                                  bool suminto,
00188                                  int vectorIndex)
00189 {
00190   int offset=0;
00191   for(int i=0; i<numIDs; ++i) {
00192     int numValues = numValuesPerID[i];
00193     if (Map().MyGID(GIDs[i])) {
00194       if (suminto) {
00195   for(int j=0; j<numValues; ++j) {
00196     SumIntoGlobalValue(GIDs[i], j, vectorIndex, values[offset+j]);
00197   }
00198       }
00199       else {
00200   for(int j=0; j<numValues; ++j) {
00201     ReplaceGlobalValue(GIDs[i], j, vectorIndex, values[offset+j]);
00202   }
00203       }
00204     }
00205     else {
00206       if (!ignoreNonLocalEntries_) {
00207   EPETRA_CHK_ERR( inputNonlocalValues(GIDs[i], numValues,
00208               &(values[offset]), suminto,
00209                                             vectorIndex) );
00210       }
00211     }
00212     offset += numValues;
00213   }
00214 
00215   return(0);
00216 }
00217 
00218 //----------------------------------------------------------------------------
00219 int Epetra_FEVector::inputNonlocalValue(int GID, double value, bool suminto,
00220                                         int vectorIndex)
00221 {
00222   int insertPoint = -1;
00223 
00224   //find offset of GID in nonlocalIDs_
00225   int offset = Epetra_Util_binary_search(GID, nonlocalIDs_, numNonlocalIDs_,
00226            insertPoint);
00227   int elemSize = Map().MaxElementSize();
00228   if (offset >= 0) {
00229     //if offset >= 0
00230     //  put value in nonlocalCoefs_[vectorIndex][offset*elemSize]
00231 
00232     offset = offset*elemSize;
00233     if (suminto) {
00234       nonlocalCoefs_[vectorIndex][offset] += value;
00235     }
00236     else {
00237       nonlocalCoefs_[vectorIndex][offset] = value;
00238     }
00239   }
00240   else {
00241     //else
00242     //  insert GID in nonlocalIDs_
00243     //  insert 1   in nonlocalElementSize_
00244     //  insert value in nonlocalCoefs_[vectorIndex]
00245 
00246     int tmp1 = numNonlocalIDs_;
00247     int tmp2 = numNonlocalIDsAlloc_;
00248     EPETRA_CHK_ERR( Epetra_Util_insert(GID, insertPoint, nonlocalIDs_,
00249                tmp1, tmp2) );
00250     --tmp1;
00251     EPETRA_CHK_ERR( Epetra_Util_insert(1, insertPoint, nonlocalElementSize_,
00252                tmp1, numNonlocalIDsAlloc_) );
00253 
00254     numNonlocalIDs_ = tmp1;
00255 
00256     //to keep nonlocalCoefs_[i] the same length for each vector in the multi-
00257     //vector, we'll insert positions for each vector even though values are
00258     //only being set for one of them...
00259     for(int i=0; i<NumVectors(); ++i) {
00260       tmp1 = numNonlocalCoefs_;
00261       tmp2 = numNonlocalCoefsAlloc_;
00262       EPETRA_CHK_ERR( Epetra_Util_insert_empty_positions(nonlocalCoefs_[i],
00263                                        tmp1, tmp2,
00264                                        insertPoint*elemSize, elemSize));
00265       for(int ii=0; ii<elemSize; ++ii) {
00266         nonlocalCoefs_[i][insertPoint*elemSize+ii] = 0.0;
00267       }
00268     }
00269     numNonlocalCoefs_ = tmp1;
00270     numNonlocalCoefsAlloc_ = tmp2;
00271 
00272     nonlocalCoefs_[vectorIndex][insertPoint*elemSize] = value;
00273   }
00274 
00275   return(0);
00276 }
00277 
00278 //----------------------------------------------------------------------------
00279 int Epetra_FEVector::inputNonlocalValues(int GID, int numValues,
00280            const double* values, bool suminto,
00281                                          int vectorIndex)
00282 {
00283   int insertPoint = -1;
00284 
00285   //find offset of GID in nonlocalIDs_
00286   int offset = Epetra_Util_binary_search(GID, nonlocalIDs_, numNonlocalIDs_,
00287            insertPoint);
00288   int elemSize = Map().MaxElementSize();
00289   if (offset >= 0) {
00290     //if offset >= 0
00291     //  put value in nonlocalCoefs_[vectorIndex][offset*elemSize]
00292 
00293     if (numValues != nonlocalElementSize_[offset]) {
00294       cerr << "Epetra_FEVector ERROR: block-size for GID " << GID << " is "
00295      << numValues<<" which doesn't match previously set block-size of "
00296      << nonlocalElementSize_[offset] << endl;
00297       return(-1);
00298     }
00299 
00300     offset = offset*elemSize;
00301 
00302     if (suminto) {
00303       for(int j=0; j<numValues; ++j) {
00304   nonlocalCoefs_[vectorIndex][offset+j] += values[j];
00305       }
00306     }
00307     else {
00308       for(int j=0; j<numValues; ++j) {
00309   nonlocalCoefs_[vectorIndex][offset+j] = values[j];
00310       }
00311     }
00312   }
00313   else {
00314     //else
00315     //  insert GID in nonlocalIDs_
00316     //  insert numValues   in nonlocalElementSize_
00317     //  insert values in nonlocalCoefs_
00318 
00319     int tmp1 = numNonlocalIDs_;
00320     int tmp2 = numNonlocalIDsAlloc_;
00321     EPETRA_CHK_ERR( Epetra_Util_insert(GID, insertPoint, nonlocalIDs_,
00322                tmp1, tmp2) );
00323     --tmp1;
00324     EPETRA_CHK_ERR( Epetra_Util_insert(numValues, insertPoint, nonlocalElementSize_,
00325                tmp1, numNonlocalIDsAlloc_) );
00326 
00327     numNonlocalIDs_ = tmp1;
00328 
00329     //to keep nonlocalCoefs_[i] the same length for each vector in the multi-
00330     //vector, we'll insert positions for each vector even though values are
00331     //only being set for one of them...
00332     for(int i=0; i<NumVectors(); ++i) {
00333       tmp1 = numNonlocalCoefs_;
00334       tmp2 = numNonlocalCoefsAlloc_;
00335       EPETRA_CHK_ERR( Epetra_Util_insert_empty_positions(nonlocalCoefs_[i],
00336                                        tmp1, tmp2,
00337                insertPoint*elemSize, elemSize));
00338       for(int ii=0; ii<elemSize; ++ii) {
00339         nonlocalCoefs_[i][insertPoint*elemSize+ii] = 0.0;
00340       }
00341     }
00342     numNonlocalCoefs_ = tmp1;
00343     numNonlocalCoefsAlloc_ = tmp2;
00344 
00345     for(int j=0; j<numValues; ++j) {
00346       nonlocalCoefs_[vectorIndex][insertPoint*elemSize+j] = values[j];
00347     }
00348   }
00349 
00350   return(0);
00351 }
00352 
00353 //----------------------------------------------------------------------------
00354 int Epetra_FEVector::GlobalAssemble(Epetra_CombineMode mode)
00355 {
00356   //In this method we need to gather all the non-local (overlapping) data
00357   //that's been input on each processor, into the (probably) non-overlapping
00358   //distribution defined by the map that 'this' vector was constructed with.
00359 
00360   //We don't need to do anything if there's only one processor or if
00361   //ignoreNonLocalEntries_ is true.
00362   if (Map().Comm().NumProc() < 2 || ignoreNonLocalEntries_) {
00363     return(0);
00364   }
00365 
00366   //First build a map that describes the data in nonlocalIDs_/nonlocalCoefs_.
00367   //We'll use the arbitrary distribution constructor of Map.
00368 
00369   Epetra_BlockMap sourceMap(-1, numNonlocalIDs_,
00370                             nonlocalIDs_, nonlocalElementSize_,
00371           Map().IndexBase(), Map().Comm());
00372 
00373   //Now build a vector to hold our nonlocalCoefs_, and to act as the source-
00374   //vector for our import operation.
00375   Epetra_MultiVector nonlocalVector(sourceMap, NumVectors());
00376 
00377   int i,j;
00378   int elemSize = Map().MaxElementSize();
00379   for(int vi=0; vi<NumVectors(); ++vi) {
00380     for(i=0; i<numNonlocalIDs_; ++i) {
00381       for(j=0; j<nonlocalElementSize_[i]; ++j) {
00382         nonlocalVector.ReplaceGlobalValue(nonlocalIDs_[i], j, vi,
00383                                           nonlocalCoefs_[vi][i*elemSize+j]);
00384       }
00385     }
00386   }
00387 
00388   Epetra_Export exporter(sourceMap, Map());
00389 
00390   EPETRA_CHK_ERR( Export(nonlocalVector, exporter, mode) );
00391 
00392   destroyNonlocalData();
00393 
00394   return(0);
00395 }
00396 
00397 //----------------------------------------------------------------------------
00398 Epetra_FEVector& Epetra_FEVector::operator=(const Epetra_FEVector& source)
00399 {
00400   Epetra_MultiVector::Assign(source);
00401 
00402   destroyNonlocalData();
00403 
00404   delete [] nonlocalCoefs_;
00405 
00406   if (source.numNonlocalIDsAlloc_ > 0) {
00407     numNonlocalIDsAlloc_ = source.numNonlocalIDsAlloc_;
00408     numNonlocalIDs_ = source.numNonlocalIDs_;
00409     nonlocalIDs_ = new int[numNonlocalIDsAlloc_];
00410     nonlocalElementSize_ = new int[numNonlocalIDsAlloc_];
00411     for(int i=0; i<numNonlocalIDs_; ++i) {
00412       int elemSize = source.nonlocalElementSize_[i];
00413       nonlocalIDs_[i] = source.nonlocalIDs_[i];
00414       nonlocalElementSize_[i] = elemSize;
00415     }
00416   }
00417 
00418   nonlocalCoefs_ = new double*[NumVectors()];
00419   for(int i=0; i<NumVectors(); ++i) nonlocalCoefs_[i] = NULL;
00420 
00421   numNonlocalCoefs_ = source.numNonlocalCoefs_;
00422   numNonlocalCoefsAlloc_ = source.numNonlocalCoefsAlloc_;
00423 
00424   if (numNonlocalCoefsAlloc_ > 0) {
00425     for(int vi=0; vi<NumVectors(); ++vi) {
00426       nonlocalCoefs_[vi] = new double[numNonlocalCoefsAlloc_];
00427       int maxelemSize = Map().MaxElementSize();
00428       for(int i=0; i<numNonlocalIDs_; ++i) {
00429         int elemSize = source.nonlocalElementSize_[i];
00430         for(int j=0; j<elemSize; ++j) {
00431           nonlocalCoefs_[vi][i*maxelemSize+j] = source.nonlocalCoefs_[vi][i*maxelemSize+j];
00432         }
00433       }
00434     }
00435   }
00436 
00437   return(*this);
00438 }
00439 
00440 //----------------------------------------------------------------------------
00441 void Epetra_FEVector::destroyNonlocalData()
00442 {
00443   if (numNonlocalIDsAlloc_ > 0) {
00444     delete [] nonlocalIDs_;
00445     delete [] nonlocalElementSize_;
00446     nonlocalIDs_ = NULL;
00447     nonlocalElementSize_ = NULL;
00448     numNonlocalIDs_ = 0;
00449     numNonlocalIDsAlloc_ = 0;
00450   }
00451 
00452   if (nonlocalCoefs_ != NULL && numNonlocalCoefsAlloc_ > 0) {
00453     for(int i=0; i<NumVectors(); ++i) {
00454       delete [] nonlocalCoefs_[i];
00455       nonlocalCoefs_[i] = NULL;
00456     }
00457 
00458     numNonlocalCoefs_ = 0;
00459     numNonlocalCoefsAlloc_ = 0;
00460   }
00461 
00462   return;
00463 }
00464 

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7