Epetra Package Browser (Single Doxygen Collection) Development
Epetra_FEVector.cpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2011 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 
00043 #include <Epetra_FEVector.h>
00044 
00045 #include <Epetra_LocalMap.h>
00046 #include <Epetra_Comm.h>
00047 #include <Epetra_Map.h>
00048 #include <Epetra_Import.h>
00049 #include <Epetra_Export.h>
00050 #include <Epetra_Util.h>
00051 #include <Epetra_IntSerialDenseVector.h>
00052 #include <Epetra_SerialDenseVector.h>
00053 
00054 #include <algorithm>
00055 
00056 //----------------------------------------------------------------------------
00057 Epetra_FEVector::Epetra_FEVector(const Epetra_BlockMap& map,
00058                                  int numVectors,
00059          bool ignoreNonLocalEntries)
00060   : Epetra_MultiVector(map, numVectors),
00061     myFirstID_(0),
00062     myNumIDs_(0),
00063     nonlocalIDs_(),
00064     nonlocalElementSize_(),
00065     nonlocalCoefs_(),
00066     nonlocalMap_(0),
00067     exporter_(0),
00068     nonlocalVector_(0),
00069     ignoreNonLocalEntries_(ignoreNonLocalEntries)
00070 {
00071   myFirstID_ = map.MinMyGID();
00072   myNumIDs_ = map.NumMyElements();
00073   nonlocalCoefs_.resize(numVectors);
00074 }
00075 
00076 //----------------------------------------------------------------------------
00077 Epetra_FEVector::Epetra_FEVector(Epetra_DataAccess CV, const Epetra_BlockMap& Map, 
00078                                  double *A, int MyLDA, int NumVectors,
00079                                  bool ignoreNonLocalEntries)
00080  : Epetra_MultiVector(CV, Map, A, MyLDA, NumVectors),
00081     myFirstID_(0),
00082     myNumIDs_(0),
00083     nonlocalIDs_(),
00084     nonlocalElementSize_(),
00085     nonlocalCoefs_(),
00086     nonlocalMap_(0),
00087     exporter_(0),
00088     nonlocalVector_(0),
00089     ignoreNonLocalEntries_(ignoreNonLocalEntries)
00090 {
00091   myFirstID_ = Map.MinMyGID();
00092   myNumIDs_ = Map.NumMyElements();
00093   nonlocalCoefs_.resize(NumVectors);
00094 }
00095 
00096 //----------------------------------------------------------------------------
00097 Epetra_FEVector::Epetra_FEVector(Epetra_DataAccess CV, const Epetra_BlockMap& Map, 
00098                                  double **ArrayOfPointers, int NumVectors,
00099                                  bool ignoreNonLocalEntries)
00100  : Epetra_MultiVector(CV, Map, ArrayOfPointers, NumVectors),
00101     myFirstID_(0),
00102     myNumIDs_(0),
00103     nonlocalIDs_(),
00104     nonlocalElementSize_(),
00105     nonlocalCoefs_(),
00106     nonlocalMap_(0),
00107     exporter_(0),
00108     nonlocalVector_(0),
00109     ignoreNonLocalEntries_(ignoreNonLocalEntries)
00110 {
00111   myFirstID_ = Map.MinMyGID();
00112   myNumIDs_ = Map.NumMyElements();
00113   nonlocalCoefs_.resize(NumVectors);
00114 }
00115 
00116 //----------------------------------------------------------------------------
00117 Epetra_FEVector::Epetra_FEVector(const Epetra_FEVector& source)
00118   : Epetra_MultiVector(source),
00119     myFirstID_(0),
00120     myNumIDs_(0),
00121     nonlocalIDs_(),
00122     nonlocalElementSize_(),
00123     nonlocalCoefs_(),
00124     nonlocalMap_(0),
00125     exporter_(0),
00126     nonlocalVector_(0),
00127     ignoreNonLocalEntries_(source.ignoreNonLocalEntries_)
00128 {
00129   *this = source;
00130 }
00131 
00132 //----------------------------------------------------------------------------
00133 Epetra_FEVector::~Epetra_FEVector()
00134 {
00135   destroyNonlocalData();
00136   destroyNonlocalMapAndExporter();
00137 
00138   nonlocalCoefs_.clear();
00139 }
00140 
00141 //----------------------------------------------------------------------------
00142 int Epetra_FEVector::SumIntoGlobalValues(int numIDs, const int* GIDs,
00143                        const double* values,
00144                                          int vectorIndex)
00145 {
00146   return( inputValues( numIDs, GIDs, values, true, vectorIndex) );
00147 }
00148 
00149 //----------------------------------------------------------------------------
00150 int Epetra_FEVector::SumIntoGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
00151                        const Epetra_SerialDenseVector& values,
00152                                          int vectorIndex)
00153 {
00154   if (GIDs.Length() != values.Length()) {
00155     return(-1);
00156   }
00157 
00158   return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), true,
00159                        vectorIndex ) );
00160 }
00161 
00162 //----------------------------------------------------------------------------
00163 int Epetra_FEVector::SumIntoGlobalValues(int numIDs, const int* GIDs,
00164            const int* numValuesPerID,
00165                        const double* values,
00166                                          int vectorIndex)
00167 {
00168   return( inputValues( numIDs, GIDs, numValuesPerID, values, true,
00169                        vectorIndex) );
00170 }
00171 
00172 //----------------------------------------------------------------------------
00173 int Epetra_FEVector::ReplaceGlobalValues(int numIDs, const int* GIDs,
00174                        const double* values,
00175                                          int vectorIndex)
00176 {
00177   return( inputValues( numIDs, GIDs, values, false,
00178                        vectorIndex) );
00179 }
00180 
00181 //----------------------------------------------------------------------------
00182 int Epetra_FEVector::ReplaceGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
00183                        const Epetra_SerialDenseVector& values,
00184                                          int vectorIndex)
00185 {
00186   if (GIDs.Length() != values.Length()) {
00187     return(-1);
00188   }
00189 
00190   return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), false,
00191                        vectorIndex) );
00192 }
00193 
00194 //----------------------------------------------------------------------------
00195 int Epetra_FEVector::inputValues(int numIDs,
00196                                  const int* GIDs,
00197                                  const double* values,
00198                                  bool suminto,
00199                                  int vectorIndex)
00200 {
00201  //Important note!! This method assumes that there is only 1 point
00202  //associated with each element (GID), and writes to offset 0 in that
00203  //GID's block.
00204 
00205   for(int i=0; i<numIDs; ++i) {
00206     if (Map().MyGID(GIDs[i])) {
00207       if (suminto) {
00208         SumIntoGlobalValue(GIDs[i], 0, vectorIndex, values[i]);
00209       }
00210       else {
00211         ReplaceGlobalValue(GIDs[i], 0, vectorIndex, values[i]);
00212       }
00213     }
00214     else {
00215       if (!ignoreNonLocalEntries_) {
00216         EPETRA_CHK_ERR( inputNonlocalValue(GIDs[i], values[i], suminto,
00217               vectorIndex) );
00218       }
00219     }
00220   }
00221 
00222   return(0);
00223 }
00224 
00225 //----------------------------------------------------------------------------
00226 int Epetra_FEVector::ReplaceGlobalValues(int numIDs, const int* GIDs,
00227            const int* numValuesPerID,
00228                        const double* values,
00229                                          int vectorIndex)
00230 {
00231   return( inputValues( numIDs, GIDs, numValuesPerID, values, false,
00232                        vectorIndex) );
00233 }
00234 
00235 //----------------------------------------------------------------------------
00236 int Epetra_FEVector::inputValues(int numIDs,
00237                                  const int* GIDs,
00238                                  const int* numValuesPerID,
00239                                  const double* values,
00240                                  bool suminto,
00241                                  int vectorIndex)
00242 {
00243   int offset=0;
00244   for(int i=0; i<numIDs; ++i) {
00245     int numValues = numValuesPerID[i];
00246     if (Map().MyGID(GIDs[i])) {
00247       if (suminto) {
00248         for(int j=0; j<numValues; ++j) {
00249           SumIntoGlobalValue(GIDs[i], j, vectorIndex, values[offset+j]);
00250         }
00251       }
00252       else {
00253         for(int j=0; j<numValues; ++j) {
00254           ReplaceGlobalValue(GIDs[i], j, vectorIndex, values[offset+j]);
00255         }
00256       }
00257     }
00258     else {
00259       if (!ignoreNonLocalEntries_) {
00260         EPETRA_CHK_ERR( inputNonlocalValues(GIDs[i], numValues,
00261               &(values[offset]), suminto,
00262               vectorIndex) );
00263       }
00264     }
00265     offset += numValues;
00266   }
00267 
00268   return(0);
00269 }
00270 
00271 //----------------------------------------------------------------------------
00272 int Epetra_FEVector::inputNonlocalValue(int GID, double value, bool suminto,
00273                                         int vectorIndex)
00274 {
00275   return inputNonlocalValues(GID, 1, &value, suminto, vectorIndex);
00276 }
00277 
00278 //----------------------------------------------------------------------------
00279 int Epetra_FEVector::inputNonlocalValues(int GID, int numValues,
00280                                          const double* values, bool suminto,
00281                                          int vectorIndex)
00282 {
00283   //find offset of GID in nonlocalIDs_
00284   std::vector<int>::iterator it = std::lower_bound(nonlocalIDs_.begin(), nonlocalIDs_.end(), GID);
00285   int offset = it - nonlocalIDs_.begin();
00286   int insertPoint = offset;
00287   if (it == nonlocalIDs_.end() || *it != GID) {
00288     offset = -1;
00289   }
00290 
00291   int elemSize = Map().MaxElementSize();
00292   if (offset >= 0) {
00293     //if offset >= 0 (meaning GID was found)
00294     //  put value in nonlocalCoefs_[vectorIndex][offset*elemSize]
00295 
00296     if (numValues != nonlocalElementSize_[offset]) {
00297       cerr << "Epetra_FEVector ERROR: block-size for GID " << GID << " is "
00298      << numValues<<" which doesn't match previously set block-size of "
00299      << nonlocalElementSize_[offset] << endl;
00300       return(-1);
00301     }
00302 
00303     offset = offset*elemSize;
00304 
00305     if (suminto) {
00306       for(int j=0; j<numValues; ++j) {
00307         nonlocalCoefs_[vectorIndex][offset+j] += values[j];
00308       }
00309     }
00310     else {
00311       for(int j=0; j<numValues; ++j) {
00312         nonlocalCoefs_[vectorIndex][offset+j] = values[j];
00313       }
00314     }
00315   }
00316   else {
00317     //else
00318     //  insert GID in nonlocalIDs_
00319     //  insert numValues   in nonlocalElementSize_
00320     //  insert values in nonlocalCoefs_
00321 
00322     nonlocalIDs_.insert(it, GID);
00323     nonlocalElementSize_.insert(nonlocalElementSize_.begin()+insertPoint, numValues);
00324 
00325     //to keep nonlocalCoefs_[i] the same length for each vector in the multi-
00326     //vector, we'll insert positions for each vector even though values are
00327     //only being set for one of them...
00328     for(int i=0; i<NumVectors(); ++i) {
00329       for(int ii=0; ii<elemSize; ++ii) {
00330         nonlocalCoefs_[i].insert(nonlocalCoefs_[i].begin()+insertPoint*elemSize+ii, 0.0);
00331       }
00332     }
00333 
00334     for(int j=0; j<numValues; ++j) {
00335       nonlocalCoefs_[vectorIndex][insertPoint*elemSize+j] = values[j];
00336     }
00337   }
00338 
00339   return(0);
00340 }
00341 
00342 //----------------------------------------------------------------------------
00343 int Epetra_FEVector::GlobalAssemble(Epetra_CombineMode mode,
00344                                     bool reuse_map_and_exporter)
00345 {
00346   //In this method we need to gather all the non-local (overlapping) data
00347   //that's been input on each processor, into the (probably) non-overlapping
00348   //distribution defined by the map that 'this' vector was constructed with.
00349 
00350   //We don't need to do anything if there's only one processor or if
00351   //ignoreNonLocalEntries_ is true.
00352   if (Map().Comm().NumProc() < 2 || ignoreNonLocalEntries_) {
00353     return(0);
00354   }
00355 
00356   if (nonlocalMap_ == 0 || !reuse_map_and_exporter) {
00357     createNonlocalMapAndExporter();
00358   }
00359 
00360   Epetra_MultiVector& nonlocalVector = *nonlocalVector_;
00361   nonlocalVector.PutScalar(0.0);
00362 
00363   int elemSize = Map().MaxElementSize();
00364   for(int vi=0; vi<NumVectors(); ++vi) {
00365     for(size_t i=0; i<nonlocalIDs_.size(); ++i) {
00366       for(int j=0; j<nonlocalElementSize_[i]; ++j) {
00367         nonlocalVector.ReplaceGlobalValue(nonlocalIDs_[i], j, vi,
00368                                           nonlocalCoefs_[vi][i*elemSize+j]);
00369       }
00370     }
00371   }
00372 
00373   EPETRA_CHK_ERR( Export(nonlocalVector, *exporter_, mode) );
00374 
00375   if (reuse_map_and_exporter) {
00376     zeroNonlocalData();
00377   }
00378   else {
00379     destroyNonlocalData();
00380   }
00381 
00382   return(0);
00383 }
00384 
00385 //----------------------------------------------------------------------------
00386 void Epetra_FEVector::createNonlocalMapAndExporter()
00387 {
00388   delete nonlocalMap_;
00389   int* nlIDptr = nonlocalIDs_.size()>0 ? &nonlocalIDs_[0] : NULL;
00390   int* nlElSzptr = nonlocalElementSize_.size()>0 ? &nonlocalElementSize_[0] : NULL;
00391   nonlocalMap_ = new Epetra_BlockMap (-1, nonlocalIDs_.size(), nlIDptr, 
00392               nlElSzptr, Map().IndexBase(), Map().Comm());
00393   delete exporter_;
00394   exporter_ = new Epetra_Export (*nonlocalMap_, Map());
00395 
00396   delete nonlocalVector_;
00397   nonlocalVector_ = new Epetra_MultiVector (*nonlocalMap_, NumVectors());
00398 }
00399 
00400 //----------------------------------------------------------------------------
00401 void Epetra_FEVector::destroyNonlocalMapAndExporter()
00402 {
00403   delete nonlocalMap_;    nonlocalMap_ = 0;
00404   delete exporter_;       exporter_ = 0;
00405   delete nonlocalVector_; nonlocalVector_ = 0;
00406 }
00407 
00408 //----------------------------------------------------------------------------
00409 Epetra_FEVector& Epetra_FEVector::operator=(const Epetra_FEVector& source)
00410 {
00411   if (this == &source) {
00412     // Don't allow self-assignment, since the allocations and
00413     // deallocations in the code below assume that source is a
00414     // different object than *this.
00415     return *this; 
00416   }
00417   // This redundantly checks for self-assignment, but the check is
00418   // inexpensive (just a pointer comparison).
00419   Epetra_MultiVector::Assign(source);
00420 
00421   nonlocalIDs_ = source.nonlocalIDs_;
00422   nonlocalElementSize_ = source.nonlocalElementSize_;
00423   nonlocalCoefs_ = source.nonlocalCoefs_;
00424 
00425   return(*this);
00426 }
00427 
00428 //----------------------------------------------------------------------------
00429 void Epetra_FEVector::zeroNonlocalData()
00430 {
00431   if (nonlocalIDs_.size() > 0) {
00432     int maxelemSize = Map().MaxElementSize();
00433     for(int vi=0; vi<NumVectors(); ++vi) {
00434       for(size_t i=0; i<nonlocalIDs_.size(); ++i) {
00435         int elemSize = nonlocalElementSize_[i];
00436         for(int j=0; j<elemSize; ++j) {
00437           nonlocalCoefs_[vi][i*maxelemSize+j] = 0.0;
00438         }
00439       }
00440     }
00441   }
00442 }
00443 
00444 //----------------------------------------------------------------------------
00445 void Epetra_FEVector::destroyNonlocalData()
00446 {
00447   nonlocalIDs_.clear();
00448   nonlocalElementSize_.clear();
00449 
00450   if (nonlocalCoefs_.size() > 0) {
00451     for(int i=0; i<NumVectors(); ++i) {
00452       nonlocalCoefs_[i].clear();
00453     }
00454   }
00455 }
00456 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines