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 2001 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 //----------------------------------------------------------------------------
00055 Epetra_FEVector::Epetra_FEVector(const Epetra_BlockMap& Map,
00056                                  int numVectors,
00057          bool ignoreNonLocalEntries)
00058   : Epetra_MultiVector(Map, numVectors),
00059     myFirstID_(0),
00060     myNumIDs_(0),
00061     nonlocalIDs_(NULL),
00062     nonlocalElementSize_(NULL),
00063     numNonlocalIDs_(0),
00064     numNonlocalIDsAlloc_(0),
00065     nonlocalCoefs_(NULL),
00066     numNonlocalCoefs_(0),
00067     numNonlocalCoefsAlloc_(0),
00068     nonlocalMap_(NULL),
00069     exporter_(NULL),
00070     nonlocalVector_(NULL),
00071     ignoreNonLocalEntries_(ignoreNonLocalEntries)
00072 {
00073   myFirstID_ = Map.MinMyGID();
00074   myNumIDs_ = Map.NumMyElements();
00075   nonlocalCoefs_ = new double*[numVectors];
00076   for(int i=0; i<numVectors; ++i) nonlocalCoefs_[i] = NULL;
00077 }
00078 
00079 //----------------------------------------------------------------------------
00080 Epetra_FEVector::Epetra_FEVector(const Epetra_FEVector& source)
00081   : Epetra_MultiVector(source),
00082     myFirstID_(0),
00083     myNumIDs_(0),
00084     nonlocalIDs_(NULL),
00085     nonlocalElementSize_(NULL),
00086     numNonlocalIDs_(0),
00087     numNonlocalIDsAlloc_(0),
00088     nonlocalCoefs_(NULL),
00089     numNonlocalCoefsAlloc_(0),
00090     nonlocalMap_(NULL),
00091     exporter_(NULL),
00092     nonlocalVector_(NULL),
00093     ignoreNonLocalEntries_(source.ignoreNonLocalEntries_)
00094 {
00095   *this = source;
00096 }
00097 
00098 //----------------------------------------------------------------------------
00099 Epetra_FEVector::~Epetra_FEVector()
00100 {
00101   destroyNonlocalData();
00102   destroyNonlocalMapAndExporter();
00103 
00104   delete [] nonlocalCoefs_;
00105   nonlocalCoefs_ = NULL;
00106 }
00107 
00108 //----------------------------------------------------------------------------
00109 int Epetra_FEVector::SumIntoGlobalValues(int numIDs, const int* GIDs,
00110                        const double* values,
00111                                          int vectorIndex)
00112 {
00113   return( inputValues( numIDs, GIDs, values, true, vectorIndex) );
00114 }
00115 
00116 //----------------------------------------------------------------------------
00117 int Epetra_FEVector::SumIntoGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
00118                        const Epetra_SerialDenseVector& values,
00119                                          int vectorIndex)
00120 {
00121   if (GIDs.Length() != values.Length()) {
00122     return(-1);
00123   }
00124 
00125   return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), true,
00126                        vectorIndex ) );
00127 }
00128 
00129 //----------------------------------------------------------------------------
00130 int Epetra_FEVector::SumIntoGlobalValues(int numIDs, const int* GIDs,
00131            const int* numValuesPerID,
00132                        const double* values,
00133                                          int vectorIndex)
00134 {
00135   return( inputValues( numIDs, GIDs, numValuesPerID, values, true,
00136                        vectorIndex) );
00137 }
00138 
00139 //----------------------------------------------------------------------------
00140 int Epetra_FEVector::ReplaceGlobalValues(int numIDs, const int* GIDs,
00141                        const double* values,
00142                                          int vectorIndex)
00143 {
00144   return( inputValues( numIDs, GIDs, values, false,
00145                        vectorIndex) );
00146 }
00147 
00148 //----------------------------------------------------------------------------
00149 int Epetra_FEVector::ReplaceGlobalValues(const Epetra_IntSerialDenseVector& GIDs,
00150                        const Epetra_SerialDenseVector& values,
00151                                          int vectorIndex)
00152 {
00153   if (GIDs.Length() != values.Length()) {
00154     return(-1);
00155   }
00156 
00157   return( inputValues( GIDs.Length(), GIDs.Values(), values.Values(), false,
00158                        vectorIndex) );
00159 }
00160 
00161 //----------------------------------------------------------------------------
00162 int Epetra_FEVector::inputValues(int numIDs,
00163                                  const int* GIDs,
00164                                  const double* values,
00165                                  bool suminto,
00166                                  int vectorIndex)
00167 {
00168  //Important note!! This method assumes that there is only 1 point
00169  //associated with each element (GID), and writes to offset 0 in that
00170  //GID's block.
00171 
00172   for(int i=0; i<numIDs; ++i) {
00173     if (Map().MyGID(GIDs[i])) {
00174       if (suminto) {
00175         SumIntoGlobalValue(GIDs[i], 0, vectorIndex, values[i]);
00176       }
00177       else {
00178         ReplaceGlobalValue(GIDs[i], 0, vectorIndex, values[i]);
00179       }
00180     }
00181     else {
00182       if (!ignoreNonLocalEntries_) {
00183         EPETRA_CHK_ERR( inputNonlocalValue(GIDs[i], values[i], suminto,
00184               vectorIndex) );
00185       }
00186     }
00187   }
00188 
00189   return(0);
00190 }
00191 
00192 //----------------------------------------------------------------------------
00193 int Epetra_FEVector::ReplaceGlobalValues(int numIDs, const int* GIDs,
00194            const int* numValuesPerID,
00195                        const double* values,
00196                                          int vectorIndex)
00197 {
00198   return( inputValues( numIDs, GIDs, numValuesPerID, values, false,
00199                        vectorIndex) );
00200 }
00201 
00202 //----------------------------------------------------------------------------
00203 int Epetra_FEVector::inputValues(int numIDs,
00204                                  const int* GIDs,
00205          const int* numValuesPerID,
00206                                  const double* values,
00207                                  bool suminto,
00208                                  int vectorIndex)
00209 {
00210   int offset=0;
00211   for(int i=0; i<numIDs; ++i) {
00212     int numValues = numValuesPerID[i];
00213     if (Map().MyGID(GIDs[i])) {
00214       if (suminto) {
00215   for(int j=0; j<numValues; ++j) {
00216     SumIntoGlobalValue(GIDs[i], j, vectorIndex, values[offset+j]);
00217   }
00218       }
00219       else {
00220   for(int j=0; j<numValues; ++j) {
00221     ReplaceGlobalValue(GIDs[i], j, vectorIndex, values[offset+j]);
00222   }
00223       }
00224     }
00225     else {
00226       if (!ignoreNonLocalEntries_) {
00227   EPETRA_CHK_ERR( inputNonlocalValues(GIDs[i], numValues,
00228               &(values[offset]), suminto,
00229                                             vectorIndex) );
00230       }
00231     }
00232     offset += numValues;
00233   }
00234 
00235   return(0);
00236 }
00237 
00238 //----------------------------------------------------------------------------
00239 int Epetra_FEVector::inputNonlocalValue(int GID, double value, bool suminto,
00240                                         int vectorIndex)
00241 {
00242   int insertPoint = -1;
00243 
00244   //find offset of GID in nonlocalIDs_
00245   int offset = Epetra_Util_binary_search(GID, nonlocalIDs_, numNonlocalIDs_,
00246            insertPoint);
00247   int elemSize = Map().MaxElementSize();
00248   if (offset >= 0) {
00249     //if offset >= 0
00250     //  put value in nonlocalCoefs_[vectorIndex][offset*elemSize]
00251 
00252     offset = offset*elemSize;
00253     if (suminto) {
00254       nonlocalCoefs_[vectorIndex][offset] += value;
00255     }
00256     else {
00257       nonlocalCoefs_[vectorIndex][offset] = value;
00258     }
00259   }
00260   else {
00261     //else
00262     //  insert GID in nonlocalIDs_
00263     //  insert 1   in nonlocalElementSize_
00264     //  insert value in nonlocalCoefs_[vectorIndex]
00265 
00266     int tmp1 = numNonlocalIDs_;
00267     int tmp2 = numNonlocalIDsAlloc_;
00268     EPETRA_CHK_ERR( Epetra_Util_insert(GID, insertPoint, nonlocalIDs_,
00269                tmp1, tmp2) );
00270     --tmp1;
00271     EPETRA_CHK_ERR( Epetra_Util_insert(1, insertPoint, nonlocalElementSize_,
00272                tmp1, numNonlocalIDsAlloc_) );
00273 
00274     numNonlocalIDs_ = tmp1;
00275 
00276     //to keep nonlocalCoefs_[i] the same length for each vector in the multi-
00277     //vector, we'll insert positions for each vector even though values are
00278     //only being set for one of them...
00279     for(int i=0; i<NumVectors(); ++i) {
00280       tmp1 = numNonlocalCoefs_;
00281       tmp2 = numNonlocalCoefsAlloc_;
00282       EPETRA_CHK_ERR( Epetra_Util_insert_empty_positions(nonlocalCoefs_[i],
00283                                        tmp1, tmp2,
00284                                        insertPoint*elemSize, elemSize));
00285       for(int ii=0; ii<elemSize; ++ii) {
00286         nonlocalCoefs_[i][insertPoint*elemSize+ii] = 0.0;
00287       }
00288     }
00289     numNonlocalCoefs_ = tmp1;
00290     numNonlocalCoefsAlloc_ = tmp2;
00291 
00292     nonlocalCoefs_[vectorIndex][insertPoint*elemSize] = value;
00293   }
00294 
00295   return(0);
00296 }
00297 
00298 //----------------------------------------------------------------------------
00299 int Epetra_FEVector::inputNonlocalValues(int GID, int numValues,
00300            const double* values, bool suminto,
00301                                          int vectorIndex)
00302 {
00303   int insertPoint = -1;
00304 
00305   //find offset of GID in nonlocalIDs_
00306   int offset = Epetra_Util_binary_search(GID, nonlocalIDs_, numNonlocalIDs_,
00307            insertPoint);
00308   int elemSize = Map().MaxElementSize();
00309   if (offset >= 0) {
00310     //if offset >= 0
00311     //  put value in nonlocalCoefs_[vectorIndex][offset*elemSize]
00312 
00313     if (numValues != nonlocalElementSize_[offset]) {
00314       cerr << "Epetra_FEVector ERROR: block-size for GID " << GID << " is "
00315      << numValues<<" which doesn't match previously set block-size of "
00316      << nonlocalElementSize_[offset] << endl;
00317       return(-1);
00318     }
00319 
00320     offset = offset*elemSize;
00321 
00322     if (suminto) {
00323       for(int j=0; j<numValues; ++j) {
00324   nonlocalCoefs_[vectorIndex][offset+j] += values[j];
00325       }
00326     }
00327     else {
00328       for(int j=0; j<numValues; ++j) {
00329   nonlocalCoefs_[vectorIndex][offset+j] = values[j];
00330       }
00331     }
00332   }
00333   else {
00334     //else
00335     //  insert GID in nonlocalIDs_
00336     //  insert numValues   in nonlocalElementSize_
00337     //  insert values in nonlocalCoefs_
00338 
00339     int tmp1 = numNonlocalIDs_;
00340     int tmp2 = numNonlocalIDsAlloc_;
00341     EPETRA_CHK_ERR( Epetra_Util_insert(GID, insertPoint, nonlocalIDs_,
00342                tmp1, tmp2) );
00343     --tmp1;
00344     EPETRA_CHK_ERR( Epetra_Util_insert(numValues, insertPoint, nonlocalElementSize_,
00345                tmp1, numNonlocalIDsAlloc_) );
00346 
00347     numNonlocalIDs_ = tmp1;
00348 
00349     //to keep nonlocalCoefs_[i] the same length for each vector in the multi-
00350     //vector, we'll insert positions for each vector even though values are
00351     //only being set for one of them...
00352     for(int i=0; i<NumVectors(); ++i) {
00353       tmp1 = numNonlocalCoefs_;
00354       tmp2 = numNonlocalCoefsAlloc_;
00355       EPETRA_CHK_ERR( Epetra_Util_insert_empty_positions(nonlocalCoefs_[i],
00356                                        tmp1, tmp2,
00357                insertPoint*elemSize, elemSize));
00358       for(int ii=0; ii<elemSize; ++ii) {
00359         nonlocalCoefs_[i][insertPoint*elemSize+ii] = 0.0;
00360       }
00361     }
00362     numNonlocalCoefs_ = tmp1;
00363     numNonlocalCoefsAlloc_ = tmp2;
00364 
00365     for(int j=0; j<numValues; ++j) {
00366       nonlocalCoefs_[vectorIndex][insertPoint*elemSize+j] = values[j];
00367     }
00368   }
00369 
00370   return(0);
00371 }
00372 
00373 //----------------------------------------------------------------------------
00374 int Epetra_FEVector::GlobalAssemble(Epetra_CombineMode mode,
00375                                     bool reuse_map_and_exporter)
00376 {
00377   //In this method we need to gather all the non-local (overlapping) data
00378   //that's been input on each processor, into the (probably) non-overlapping
00379   //distribution defined by the map that 'this' vector was constructed with.
00380 
00381   //We don't need to do anything if there's only one processor or if
00382   //ignoreNonLocalEntries_ is true.
00383   if (Map().Comm().NumProc() < 2 || ignoreNonLocalEntries_) {
00384     return(0);
00385   }
00386 
00387   if (nonlocalMap_ == NULL || !reuse_map_and_exporter) {
00388     createNonlocalMapAndExporter();
00389   }
00390 
00391   Epetra_MultiVector& nonlocalVector = *nonlocalVector_;
00392   nonlocalVector.PutScalar(0.0);
00393 
00394   int elemSize = Map().MaxElementSize();
00395   for(int vi=0; vi<NumVectors(); ++vi) {
00396     for(int i=0; i<numNonlocalIDs_; ++i) {
00397       for(int j=0; j<nonlocalElementSize_[i]; ++j) {
00398         nonlocalVector.ReplaceGlobalValue(nonlocalIDs_[i], j, vi,
00399                                           nonlocalCoefs_[vi][i*elemSize+j]);
00400       }
00401     }
00402   }
00403 
00404   EPETRA_CHK_ERR( Export(nonlocalVector, *exporter_, mode) );
00405 
00406   zeroNonlocalData();
00407 
00408   return(0);
00409 }
00410 
00411 //----------------------------------------------------------------------------
00412 void Epetra_FEVector::createNonlocalMapAndExporter()
00413 {
00414   delete nonlocalMap_;
00415   nonlocalMap_ = new Epetra_BlockMap(-1, numNonlocalIDs_, nonlocalIDs_,
00416                                nonlocalElementSize_, Map().IndexBase(), Map().Comm());
00417   delete exporter_;
00418   exporter_ = new Epetra_Export(*nonlocalMap_, Map());
00419 
00420   delete nonlocalVector_;
00421   nonlocalVector_ = new Epetra_MultiVector(*nonlocalMap_, NumVectors());
00422 }
00423 
00424 //----------------------------------------------------------------------------
00425 void Epetra_FEVector::destroyNonlocalMapAndExporter()
00426 {
00427   delete nonlocalMap_;
00428   delete exporter_;
00429   delete nonlocalVector_;
00430 }
00431 
00432 //----------------------------------------------------------------------------
00433 Epetra_FEVector& Epetra_FEVector::operator=(const Epetra_FEVector& source)
00434 {
00435   Epetra_MultiVector::Assign(source);
00436 
00437   destroyNonlocalData();
00438 
00439   delete [] nonlocalCoefs_;
00440 
00441   if (source.numNonlocalIDsAlloc_ > 0) {
00442     numNonlocalIDsAlloc_ = source.numNonlocalIDsAlloc_;
00443     numNonlocalIDs_ = source.numNonlocalIDs_;
00444     nonlocalIDs_ = new int[numNonlocalIDsAlloc_];
00445     nonlocalElementSize_ = new int[numNonlocalIDsAlloc_];
00446     for(int i=0; i<numNonlocalIDs_; ++i) {
00447       int elemSize = source.nonlocalElementSize_[i];
00448       nonlocalIDs_[i] = source.nonlocalIDs_[i];
00449       nonlocalElementSize_[i] = elemSize;
00450     }
00451   }
00452 
00453   nonlocalCoefs_ = new double*[NumVectors()];
00454   for(int i=0; i<NumVectors(); ++i) nonlocalCoefs_[i] = NULL;
00455 
00456   numNonlocalCoefs_ = source.numNonlocalCoefs_;
00457   numNonlocalCoefsAlloc_ = source.numNonlocalCoefsAlloc_;
00458 
00459   if (numNonlocalCoefsAlloc_ > 0) {
00460     for(int vi=0; vi<NumVectors(); ++vi) {
00461       nonlocalCoefs_[vi] = new double[numNonlocalCoefsAlloc_];
00462       int maxelemSize = Map().MaxElementSize();
00463       for(int i=0; i<numNonlocalIDs_; ++i) {
00464         int elemSize = source.nonlocalElementSize_[i];
00465         for(int j=0; j<elemSize; ++j) {
00466           nonlocalCoefs_[vi][i*maxelemSize+j] = source.nonlocalCoefs_[vi][i*maxelemSize+j];
00467         }
00468       }
00469     }
00470   }
00471 
00472   return(*this);
00473 }
00474 
00475 //----------------------------------------------------------------------------
00476 void Epetra_FEVector::zeroNonlocalData()
00477 {
00478   if (numNonlocalCoefsAlloc_ > 0) {
00479     int maxelemSize = Map().MaxElementSize();
00480     for(int vi=0; vi<NumVectors(); ++vi) {
00481       for(int i=0; i<numNonlocalIDs_; ++i) {
00482         int elemSize = nonlocalElementSize_[i];
00483         for(int j=0; j<elemSize; ++j) {
00484           nonlocalCoefs_[vi][i*maxelemSize+j] = 0.0;
00485         }
00486       }
00487     }
00488   }
00489 }
00490 
00491 //----------------------------------------------------------------------------
00492 void Epetra_FEVector::destroyNonlocalData()
00493 {
00494   if (numNonlocalIDsAlloc_ > 0) {
00495     delete [] nonlocalIDs_;
00496     delete [] nonlocalElementSize_;
00497     nonlocalIDs_ = NULL;
00498     nonlocalElementSize_ = NULL;
00499     numNonlocalIDs_ = 0;
00500     numNonlocalIDsAlloc_ = 0;
00501   }
00502 
00503   if (nonlocalCoefs_ != NULL && numNonlocalCoefsAlloc_ > 0) {
00504     for(int i=0; i<NumVectors(); ++i) {
00505       delete [] nonlocalCoefs_[i];
00506       nonlocalCoefs_[i] = NULL;
00507     }
00508 
00509     numNonlocalCoefs_ = 0;
00510     numNonlocalCoefsAlloc_ = 0;
00511   }
00512 }
00513 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines