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 //Initial implementation of Epetra_FEVector. At this point it is
00043 //capable of accepting overlapping data and then gathering that data onto
00044 //the "owning" processors when GlobalAssemble() is called.
00045 //
00046 //Issues to be resolved:
00047 //1. Assumptions are currently made implicitly, that ElementSize==1 for the
00048 //Epetra_BlockMap that's provided at construction. Code needs to be
00049 //generalized for varying element-sizes.
00050 //
00051 //2. Implementation of SumIntoGlobalValues() and ReplaceGlobalValues() is
00052 //currently inefficient (probably). Code needs to be optimized.
00053 //
00054 
00055 //----------------------------------------------------------------------------
00056 Epetra_FEVector::Epetra_FEVector(const Epetra_BlockMap& Map,
00057          bool ignoreNonLocalEntries)
00058   : Epetra_MultiVector(Map, 1),
00059     myFirstID_(0),
00060     myNumIDs_(0),
00061     myCoefs_(NULL),
00062     nonlocalIDs_(NULL),
00063     nonlocalElementSize_(NULL),
00064     numNonlocalIDs_(0),
00065     allocatedNonlocalLength_(0),
00066     nonlocalCoefs_(NULL),
00067     ignoreNonLocalEntries_(ignoreNonLocalEntries)
00068 {
00069   myFirstID_ = Map.MinMyGID();
00070   myNumIDs_ = Map.NumMyElements();
00071 
00072   //Currently we impose the restriction that NumVectors==1, so we won't
00073   //need the LDA argument when calling ExtractView. Hence the "dummy" arg.
00074   int dummy;
00075   ExtractView(&myCoefs_, &dummy);
00076 }
00077 
00078 //----------------------------------------------------------------------------
00079 Epetra_FEVector::Epetra_FEVector(const Epetra_FEVector& source)
00080   : Epetra_MultiVector(source.Map(), 1),
00081     myFirstID_(0),
00082     myNumIDs_(0),
00083     myCoefs_(NULL),
00084     nonlocalIDs_(NULL),
00085     nonlocalElementSize_(NULL),
00086     numNonlocalIDs_(0),
00087     allocatedNonlocalLength_(0),
00088     nonlocalCoefs_(NULL),
00089     ignoreNonLocalEntries_(source.ignoreNonLocalEntries_)
00090 {
00091   *this = source;
00092 }
00093 
00094 //----------------------------------------------------------------------------
00095 Epetra_FEVector::~Epetra_FEVector()
00096 {
00097   destroyNonlocalData();
00098 }
00099 
00100 //----------------------------------------------------------------------------
00101 int Epetra_FEVector::SumIntoGlobalValues(int numIDs, const int* GIDs,
00102                        const double* values)
00103 {
00104   return( inputValues( numIDs, GIDs, values, true) );
00105 }
00106 
00107 //----------------------------------------------------------------------------
00108 int Epetra_FEVector::SumIntoGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
00109                        const Epetra_SerialDenseVector& values)
00110 {
00111   if (GIDs.Length() != values.Length()) {
00112     return(-1);
00113   }
00114 
00115   return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), true) );
00116 }
00117 
00118 //----------------------------------------------------------------------------
00119 int Epetra_FEVector::SumIntoGlobalValues(int numIDs, const int* GIDs,
00120            const int* numValuesPerID,
00121                        const double* values)
00122 {
00123   return( inputValues( numIDs, GIDs, numValuesPerID, values, true) );
00124 }
00125 
00126 //----------------------------------------------------------------------------
00127 int Epetra_FEVector::ReplaceGlobalValues(int numIDs, const int* GIDs,
00128                        const double* values)
00129 {
00130   return( inputValues( numIDs, GIDs, values, false) );
00131 }
00132 
00133 //----------------------------------------------------------------------------
00134 int Epetra_FEVector::ReplaceGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
00135                        const Epetra_SerialDenseVector& values)
00136 {
00137   if (GIDs.Length() != values.Length()) {
00138     return(-1);
00139   }
00140 
00141   return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), false) );
00142 }
00143 
00144 //----------------------------------------------------------------------------
00145 int Epetra_FEVector::ReplaceGlobalValues(int numIDs, const int* GIDs,
00146            const int* numValuesPerID,
00147                        const double* values)
00148 {
00149   return( inputValues( numIDs, GIDs, numValuesPerID, values, false) );
00150 }
00151 
00152 //----------------------------------------------------------------------------
00153 int Epetra_FEVector::inputValues(int numIDs,
00154                                  const int* GIDs,
00155                                  const double* values,
00156                                  bool accumulate)
00157 {
00158  //Important note!! This method assumes that there is only 1 point
00159  //associated with each element.
00160 
00161   for(int i=0; i<numIDs; ++i) {
00162     if (Map().MyGID(GIDs[i])) {
00163       if (accumulate) {
00164         SumIntoGlobalValue(GIDs[i], 0, 0, values[i]);
00165       }
00166       else {
00167         ReplaceGlobalValue(GIDs[i], 0, 0, values[i]);
00168       }
00169     }
00170     else {
00171       if (!ignoreNonLocalEntries_) {
00172   EPETRA_CHK_ERR( inputNonlocalValue(GIDs[i], values[i], accumulate) );
00173       }
00174     }
00175   }
00176 
00177   return(0);
00178 }
00179 
00180 //----------------------------------------------------------------------------
00181 int Epetra_FEVector::inputValues(int numIDs,
00182                                  const int* GIDs,
00183          const int* numValuesPerID,
00184                                  const double* values,
00185                                  bool accumulate)
00186 {
00187   int offset=0;
00188   for(int i=0; i<numIDs; ++i) {
00189     int numValues = numValuesPerID[i];
00190     if (Map().MyGID(GIDs[i])) {
00191       if (accumulate) {
00192   for(int j=0; j<numValues; ++j) {
00193     SumIntoGlobalValue(GIDs[i], j, 0, values[offset+j]);
00194   }
00195       }
00196       else {
00197   for(int j=0; j<numValues; ++j) {
00198     ReplaceGlobalValue(GIDs[i], j, 0, values[offset+j]);
00199   }
00200       }
00201     }
00202     else {
00203       if (!ignoreNonLocalEntries_) {
00204   EPETRA_CHK_ERR( inputNonlocalValues(GIDs[i], numValues,
00205               &(values[offset]), accumulate) );
00206       }
00207     }
00208     offset += numValues;
00209   }
00210 
00211   return(0);
00212 }
00213 
00214 //----------------------------------------------------------------------------
00215 int Epetra_FEVector::inputNonlocalValue(int GID, double value, bool accumulate)
00216 {
00217   int insertPoint = -1;
00218 
00219   //find offset of GID in nonlocalIDs_
00220   int offset = Epetra_Util_binary_search(GID, nonlocalIDs_, numNonlocalIDs_,
00221            insertPoint);
00222   if (offset >= 0) {
00223     //if offset >= 0
00224     //  put value in nonlocalCoefs_[offset][0]
00225 
00226     if (accumulate) {
00227       nonlocalCoefs_[offset][0] += value;
00228     }
00229     else {
00230       nonlocalCoefs_[offset][0] = value;
00231     }
00232   }
00233   else {
00234     //else
00235     //  insert GID in nonlocalIDs_
00236     //  insert 1   in nonlocalElementSize_
00237     //  insert value in nonlocalCoefs_
00238 
00239     int tmp1 = numNonlocalIDs_;
00240     int tmp2 = allocatedNonlocalLength_;
00241     int tmp3 = allocatedNonlocalLength_;
00242     EPETRA_CHK_ERR( Epetra_Util_insert(GID, insertPoint, nonlocalIDs_,
00243                tmp1, tmp2) );
00244     --tmp1;
00245     EPETRA_CHK_ERR( Epetra_Util_insert(1, insertPoint, nonlocalElementSize_,
00246                tmp1, tmp3) );
00247     double* values = new double[1];
00248     values[0] = value;
00249     EPETRA_CHK_ERR( Epetra_Util_insert(values, insertPoint, nonlocalCoefs_,
00250                numNonlocalIDs_, allocatedNonlocalLength_) );
00251   }
00252 
00253   return(0);
00254 }
00255 
00256 //----------------------------------------------------------------------------
00257 int Epetra_FEVector::inputNonlocalValues(int GID, int numValues,
00258            const double* values, bool accumulate)
00259 {
00260   int insertPoint = -1;
00261 
00262   //find offset of GID in nonlocalIDs_
00263   int offset = Epetra_Util_binary_search(GID, nonlocalIDs_, numNonlocalIDs_,
00264            insertPoint);
00265   if (offset >= 0) {
00266     //if offset >= 0
00267     //  put value in nonlocalCoefs_[offset][0]
00268 
00269     if (numValues != nonlocalElementSize_[offset]) {
00270       cerr << "Epetra_FEVector ERROR: block-size for GID " << GID << " is "
00271      << numValues<<" which doesn't match previously set block-size of "
00272      << nonlocalElementSize_[offset] << endl;
00273       return(-1);
00274     }
00275 
00276     if (accumulate) {
00277       for(int j=0; j<numValues; ++j) {
00278   nonlocalCoefs_[offset][j] += values[j];
00279       }
00280     }
00281     else {
00282       for(int j=0; j<numValues; ++j) {
00283   nonlocalCoefs_[offset][j] = values[j];
00284       }
00285     }
00286   }
00287   else {
00288     //else
00289     //  insert GID in nonlocalIDs_
00290     //  insert numValues   in nonlocalElementSize_
00291     //  insert values in nonlocalCoefs_
00292 
00293     int tmp1 = numNonlocalIDs_;
00294     int tmp2 = allocatedNonlocalLength_;
00295     int tmp3 = allocatedNonlocalLength_;
00296     EPETRA_CHK_ERR( Epetra_Util_insert(GID, insertPoint, nonlocalIDs_,
00297                tmp1, tmp2) );
00298     --tmp1;
00299     EPETRA_CHK_ERR( Epetra_Util_insert(numValues, insertPoint, nonlocalElementSize_,
00300                tmp1, tmp3) );
00301     double* newvalues = new double[numValues];
00302     for(int j=0; j<numValues; ++j) {
00303       newvalues[j] = values[j];
00304     }
00305     EPETRA_CHK_ERR( Epetra_Util_insert(newvalues, insertPoint, nonlocalCoefs_,
00306                numNonlocalIDs_, allocatedNonlocalLength_) );
00307   }
00308 
00309   return(0);
00310 }
00311 
00312 //----------------------------------------------------------------------------
00313 int Epetra_FEVector::GlobalAssemble(Epetra_CombineMode mode)
00314 {
00315   //In this method we need to gather all the non-local (overlapping) data
00316   //that's been input on each processor, into the (probably) non-overlapping
00317   //distribution defined by the map that 'this' vector was constructed with.
00318 
00319   //We don't need to do anything if there's only one processor or if
00320   //ignoreNonLocalEntries_ is true.
00321   if (Map().Comm().NumProc() < 2 || ignoreNonLocalEntries_) {
00322     return(0);
00323   }
00324 
00325   //First build a map that describes the data in nonlocalIDs_/nonlocalCoefs_.
00326   //We'll use the arbitrary distribution constructor of Map.
00327 
00328   Epetra_BlockMap sourceMap(-1, numNonlocalIDs_,
00329                             nonlocalIDs_, nonlocalElementSize_,
00330           Map().IndexBase(), Map().Comm());
00331 
00332   //Now build a vector to hold our nonlocalCoefs_, and to act as the source-
00333   //vector for our import operation.
00334   Epetra_MultiVector nonlocalVector(sourceMap, 1);
00335 
00336   int i,j;
00337   for(i=0; i<numNonlocalIDs_; ++i) {
00338     for(j=0; j<nonlocalElementSize_[i]; ++j) {
00339       nonlocalVector.ReplaceGlobalValue(nonlocalIDs_[i], j, 0,
00340           nonlocalCoefs_[i][j]);
00341     }
00342   }
00343 
00344   Epetra_Export exporter(sourceMap, Map());
00345 
00346   EPETRA_CHK_ERR( Export(nonlocalVector, exporter, mode) );
00347 
00348   destroyNonlocalData();
00349 
00350   return(0);
00351 }
00352 
00353 //----------------------------------------------------------------------------
00354 Epetra_FEVector& Epetra_FEVector::operator=(const Epetra_FEVector& source)
00355 {
00356   Epetra_MultiVector::Assign(source);
00357 
00358   destroyNonlocalData();
00359 
00360   if (source.allocatedNonlocalLength_ > 0) {
00361     allocatedNonlocalLength_ = source.allocatedNonlocalLength_;
00362     numNonlocalIDs_ = source.numNonlocalIDs_;
00363     nonlocalIDs_ = new int[allocatedNonlocalLength_];
00364     nonlocalElementSize_ = new int[allocatedNonlocalLength_];
00365     nonlocalCoefs_ = new double*[allocatedNonlocalLength_];
00366     for(int i=0; i<numNonlocalIDs_; ++i) {
00367       int elemSize = source.nonlocalElementSize_[i];
00368       nonlocalCoefs_[i] = new double[elemSize];
00369       nonlocalIDs_[i] = source.nonlocalIDs_[i];
00370       nonlocalElementSize_[i] = elemSize;
00371       for(int j=0; j<elemSize; ++j) {
00372   nonlocalCoefs_[i][j] = source.nonlocalCoefs_[i][j];
00373       }
00374     }
00375   }
00376 
00377   return(*this);
00378 }
00379 
00380 //----------------------------------------------------------------------------
00381 void Epetra_FEVector::destroyNonlocalData()
00382 {
00383   if (allocatedNonlocalLength_ > 0) {
00384     delete [] nonlocalIDs_;
00385     delete [] nonlocalElementSize_;
00386     nonlocalIDs_ = NULL;
00387     nonlocalElementSize_ = NULL;
00388     for(int i=0; i<numNonlocalIDs_; ++i) {
00389       delete [] nonlocalCoefs_[i];
00390     }
00391     delete [] nonlocalCoefs_;
00392     nonlocalCoefs_ = NULL;
00393     numNonlocalIDs_ = 0;
00394     allocatedNonlocalLength_ = 0;
00395   }
00396   return;
00397 }
00398 

Generated on Thu Sep 18 12:37:57 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1