FEI Version of the Day
FEI_tester.cpp
00001 /*
00002 // @HEADER
00003 // ************************************************************************
00004 //             FEI: Finite Element Interface to Linear Solvers
00005 //                  Copyright (2005) Sandia Corporation.
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, the
00008 // U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Alan Williams (william@sandia.gov) 
00038 //
00039 // ************************************************************************
00040 // @HEADER
00041 */
00042 
00043 
00044 #include <fei_sstream.hpp>
00045 #include <fei_fstream.hpp>
00046 
00047 #include <test_utils/fei_test_utils.hpp>
00048 
00049 #include <test_utils/FEI_tester.hpp>
00050 
00051 #include <fei_LinearSystemCore.hpp>
00052 #include <fei_LibraryWrapper.hpp>
00053 #include <snl_fei_Utils.hpp>
00054 
00055 #include <fei_FEI_Impl.hpp>
00056 
00057 #include <test_utils/LibraryFactory.hpp>
00058 
00059 #ifdef HAVE_FEI_FETI
00060 #include <FETI_DP_FiniteElementData.h>
00061 #endif
00062 
00063 #include <test_utils/DataReader.hpp>
00064 #include <test_utils/SolnCheck.hpp>
00065 
00066 #undef fei_file
00067 #define fei_file "FEI_tester.cpp"
00068 
00069 #include <fei_ErrMacros.hpp>
00070 
00071 //----------------------------------------------------------------------------
00072 FEI_tester::FEI_tester(fei::SharedPtr<DataReader> data_reader,
00073            MPI_Comm comm, int localProc, int numProcs, bool useNewFEI)
00074   : comm_(comm),
00075     fei_(),
00076     wrapper_(),
00077     data_(data_reader),
00078     localProc_(localProc),
00079     numProcs_(numProcs),
00080     numMatrices(1),
00081     matrixIDs(NULL),
00082     numRHSs(1),
00083     rhsIDs(NULL),
00084     useNewFEI_(useNewFEI)
00085 {
00086 }
00087 
00088 //----------------------------------------------------------------------------
00089 FEI_tester::~FEI_tester()
00090 {
00091   delete [] matrixIDs;
00092   delete [] rhsIDs;
00093 }
00094 
00095 //----------------------------------------------------------------------------
00096 int FEI_tester::testInitialization()
00097 {
00098   if (data_.get() == NULL) {
00099     ERReturn(-1);
00100   }
00101 
00102   CHK_ERR( createFEIinstance(data_->solverLibraryName_.c_str()) );
00103 
00104   const char* feiVersionString;
00105   CHK_ERR( fei_->version(feiVersionString) );
00106 
00107   FEI_COUT << "FEI version: " << feiVersionString << FEI_ENDL;
00108 
00109   fei_->parameters(data_->numParams_, data_->paramStrings_);
00110 
00111   //Now we check the solveType. A regular Ax=b solve corresponds to 
00112   //solveType==FEI_SINGLE_SYSTEM. The aggregate stuff (solveType==
00113   //FEI_AGGREGATE_SUM) is for power users who
00114   //want to assemble more than one matrix system and solve a linear
00115   //combination of them, an aggregate system.
00116 
00117   if (data_->solveType_ == FEI_AGGREGATE_SUM) {
00118     CHK_ERR( setIDlists());
00119   }
00120 
00121   CHK_ERR( initializationPhase() );
00122 
00123   int numBlkActNodes;
00124   for(int i=0; i<data_->numElemBlocks_; ++i) {
00125     ElemBlock& eblk = data_->elemBlocks_[i];
00126     int elemBlockID = eblk.blockID_;
00127     CHK_ERR( fei_->getNumBlockActNodes(elemBlockID, numBlkActNodes) );
00128   }
00129 
00130   //************** Initialization Phase is now complete *****************
00131 
00132   return(0);
00133 }
00134 
00135 //----------------------------------------------------------------------------
00136 int FEI_tester::testLoading()
00137 {
00138   CHK_ERR(fei_->resetRHSVector());
00139   CHK_ERR(fei_->resetMatrix());    //try out all of these reset functions,
00140   CHK_ERR(fei_->resetSystem());    //just for the sake of coverage...
00141 
00142   // ************ Load Phase (delegated to a function) ***************
00143 
00144   //obviously only one of these load-phases should be
00145   //performed.
00146 
00147   if (data_->solveType_ == FEI_SINGLE_SYSTEM) {
00148     CHK_ERR( normalLoadPhase());
00149   }
00150   if (data_->solveType_ == FEI_AGGREGATE_SUM) {
00151     CHK_ERR( aggregateLoadPhase());
00152   }
00153 
00154   CHK_ERR( fei_->loadComplete() );
00155 
00156   CHK_ERR( exerciseResidualNorm() );
00157 
00158   //**** let's try out the 'put' functions. ******************
00159   CHK_ERR( exercisePutFunctions() );
00160 
00161   return(0);
00162 }
00163 
00164 //----------------------------------------------------------------------------
00165 int FEI_tester::testSolve()
00166 {
00167   //If solverLibraryName is TEST_LSC, then we're not running a real solver so
00168   //just return.
00169   std::string sname(data_->solverLibraryName_);
00170   if (sname == "TEST_LSC") {
00171     return( 0 );
00172   }
00173 
00174   int status;
00175   int err = fei_->solve(status);
00176 
00177   //FEI_COUT << "fei->solve, err = " << err << ", status = " << status << FEI_ENDL;
00178 
00179   if (err != 0 || status != 0) {
00180     FEI_COUT << "!!!! solve returned: err: "<<err<<", status: "<<status<<FEI_ENDL;
00181     return(err);
00182   }
00183 
00184   if (localProc_ == 0) {
00185     int itersTaken;
00186     CHK_ERR( fei_->iterations(itersTaken));
00187     //FEI_COUT << "iterations: " << itersTaken << FEI_ENDL;
00188   }
00189 
00190   CHK_ERR( exerciseResidualNorm() );
00191 
00192   return(0);
00193 }
00194 
00195 //----------------------------------------------------------------------------
00196 void FEI_tester::dumpMatrixFiles()
00197 {
00198 }
00199 
00200 //----------------------------------------------------------------------------
00201 void FEI_tester::setParameter(const char*)
00202 {
00203 }
00204 
00205 //----------------------------------------------------------------------------
00206 int FEI_tester::testCheckResult()
00207 {
00208   //If solverLibraryName is TEST_LSC, then we're not running a real solver so
00209   //just check the matrix.
00210   std::string sname(data_->solverLibraryName_);
00211   if (sname == "TEST_LSC") {
00212     return( lsc_matrix_check() );
00213   }
00214 
00215   CHK_ERR( save_block_node_soln(*data_, *fei_, data_->solnFileName_.c_str(),
00216         numProcs_, localProc_, 1));
00217 
00218   CHK_ERR( save_block_elem_soln(*data_, *fei_, data_->solnFileName_.c_str(),
00219         numProcs_, localProc_, 1));
00220 
00221   CHK_ERR( save_multiplier_soln(*data_, *fei_, data_->solnFileName_.c_str(),
00222         numProcs_, localProc_, 1));
00223 
00224   int err = SolnCheck::checkSolution(localProc_, numProcs_, data_->solnFileName_.c_str(),
00225              data_->checkFileName_.c_str(), "node", 1);
00226 
00227   err += SolnCheck::checkSolution(localProc_, numProcs_, data_->solnFileName_.c_str(),
00228           data_->checkFileName_.c_str(), "elem", 1);
00229 
00230   err += SolnCheck::checkSolution(localProc_, numProcs_, data_->solnFileName_.c_str(),
00231           data_->checkFileName_.c_str(), "mult", 1);
00232   int globalErr = err;
00233 #ifndef FEI_SER
00234   if (MPI_SUCCESS != MPI_Allreduce(&err, &globalErr, 1, MPI_INT, MPI_SUM,
00235            comm_)) return(-1);
00236 #endif
00237   if (globalErr != 0) {
00238     ERReturn(-1);
00239   }
00240 
00241   return(0);
00242 }
00243 
00244 //----------------------------------------------------------------------------
00245 int FEI_tester::createFEIinstance(const char* solverName)
00246 {
00247   try {
00248     wrapper_ = fei::create_LibraryWrapper(comm_, solverName);
00249   }
00250   catch(std::runtime_error& exc) {
00251     fei::console_out() << exc.what()<<FEI_ENDL;
00252     return(1);
00253   }
00254 
00255   if (wrapper_.get() == NULL) ERReturn(-1);
00256 
00257   if (useNewFEI_) {
00258     fei_.reset(new fei::FEI_Impl(wrapper_, comm_, 0));
00259   }
00260   else {
00261     fei_.reset(new FEI_Implementation(wrapper_, comm_, 0));
00262   }
00263 
00264   if (fei_.get() == NULL) ERReturn(-1);
00265 
00266   return(0);
00267 }
00268 
00269 //----------------------------------------------------------------------------
00270 int FEI_tester::setIDlists()
00271 {
00272   snl_fei::getIntParamValue("numMatrices",
00273           data_->numParams_,
00274           data_->paramStrings_,
00275           numMatrices);
00276 
00277   matrixIDs = new int[numMatrices];
00278   numRHSs = 1;
00279   rhsIDs = new int[1];
00280   rhsIDs[0] = 0;
00281 
00282   for(int i=0; i<numMatrices; i++) {
00283     matrixIDs[i] = i;
00284   }
00285 
00286   CHK_ERR(fei_->setIDLists(numMatrices, matrixIDs, numRHSs, rhsIDs));
00287   return(0);
00288 }
00289 
00290 //----------------------------------------------------------------------------
00291 int FEI_tester::initializationPhase()
00292 {
00293   if (data_->solveType_ != FEI_SINGLE_SYSTEM &&
00294       data_->solveType_ != FEI_AGGREGATE_SUM) {
00295     FEI_COUT << "FEI_tester: bad solveType: " << data_->solveType_ << FEI_ENDL;
00296     return(-1);
00297   }
00298 
00299   CHK_ERR( fei_->setSolveType(data_->solveType_));
00300 
00301   CHK_ERR(fei_->initFields(data_->numFields_, data_->fieldSizes_, data_->fieldIDs_));
00302 
00303   int i;
00304   for(i=0; i<data_->numElemBlocks_; i++) {
00305     ElemBlock& block = data_->elemBlocks_[i];
00306 
00307     CHK_ERR(fei_->initElemBlock(block.blockID_,
00308              block.numElements_,
00309              block.numNodesPerElement_,
00310              block.numFieldsPerNode_,
00311              block.nodalFieldIDs_,
00312              block.numElemDOF_,
00313              block.elemDOFFieldIDs_,
00314              block.interleaveStrategy_) );
00315 
00316     for(int el=0; el<block.numElements_; el++) {
00317       CHK_ERR(fei_->initElem(block.blockID_, 
00318           block.elemIDs_[el],
00319           block.elemConn_[el]));
00320     }
00321   }
00322 
00323   for(i=0; i<data_->numSharedNodeSets_; i++) {
00324     CommNodeSet& shNodeSet = data_->sharedNodeSets_[i];
00325 
00326     CHK_ERR(fei_->initSharedNodes(shNodeSet.numNodes_, shNodeSet.nodeIDs_,
00327          shNodeSet.procsPerNode_, shNodeSet.procs_));
00328   }
00329 
00330   //********* Initialize any slave variables *****************************
00331 
00332   for(i=0; i<data_->numSlaveVars_; i++) {
00333     CRSet& crSet = data_->slaveVars_[i];
00334 
00335     CHK_ERR( fei_->initSlaveVariable(crSet.slaveNodeID_,
00336             crSet.slaveFieldID_,
00337             crSet.slaveOffset_,
00338             crSet.numNodes_,
00339             crSet.nodeIDs_[0],
00340             crSet.fieldIDs_,
00341             crSet.weights_,
00342             crSet.values_[0]) );
00343   }
00344 
00345   //*********** Initialize Constraint Relation Equations *****************
00346 
00347   for(i=0; i<data_->numCRMultSets_; i++) {
00348     CRSet& crSet = data_->crMultSets_[i];
00349 
00350     for(int j=0; j<1; j++) {
00351       CHK_ERR(fei_->initCRMult(crSet.numNodes_,
00352             crSet.nodeIDs_[j],
00353             crSet.fieldIDs_,
00354             crSet.crID_));
00355     }
00356   }
00357 
00358   for(i=0; i<data_->numCRPenSets_; i++) {
00359     CRSet& crSet = data_->crPenSets_[i];
00360 
00361     for(int j=0; j<1; j++) {
00362       CHK_ERR(fei_->initCRPen(crSet.numNodes_,
00363            crSet.nodeIDs_[j],
00364            crSet.fieldIDs_,
00365            crSet.crID_));
00366     }
00367   }
00368 
00369   CHK_ERR(fei_->initComplete());
00370   return(0);
00371 }
00372 
00373 //----------------------------------------------------------------------------
00374 int FEI_tester::normalLoadPhase()
00375 {
00376   int i;
00377 
00378   //*************** Boundary Condition Node Sets *************************
00379 
00380   for(i=0; i<data_->numBCNodeSets_; i++) {
00381     BCNodeSet& bcSet = data_->bcNodeSets_[i];
00382 
00383     CHK_ERR(fei_->loadNodeBCs(bcSet.numNodes_,
00384            bcSet.nodeIDs_,
00385            bcSet.fieldID_,
00386            bcSet.offsetsIntoField_,
00387            bcSet.prescribed_values_));
00388   }
00389 
00390    for(i=0; i<data_->numElemBlocks_; i++) {
00391       ElemBlock& block = data_->elemBlocks_[i];
00392 
00393       for(int el=0; el<block.numElements_; el++) {
00394 
00395          CHK_ERR(fei_->sumInElemMatrix(block.blockID_,
00396                                 block.elemIDs_[el],
00397                                 block.elemConn_[el],
00398                                 block.elemStiff_[el],
00399                                 block.elemFormat_));
00400       }
00401    }
00402 
00403    for(i=0; i<data_->numElemBlocks_; i++) {
00404       ElemBlock& block = data_->elemBlocks_[i];
00405 
00406       for(int el=0; el<block.numElements_; el++) {
00407 
00408          CHK_ERR(fei_->sumInElemRHS(block.blockID_,
00409                                 block.elemIDs_[el],
00410                                 block.elemConn_[el],
00411                                 block.elemLoad_[el]));
00412       }
00413    }
00414 
00415    //******** Load Constraint Relation Equations ***********************
00416 
00417    for(i=0; i<data_->numCRMultSets_; i++) {
00418       CRSet& crSet = data_->crMultSets_[i];
00419 
00420       for(int j=0; j<1; j++) {
00421          CHK_ERR(fei_->loadCRMult(crSet.crID_,
00422                                  crSet.numNodes_,
00423                                  crSet.nodeIDs_[j],
00424                                  crSet.fieldIDs_,
00425                                  crSet.weights_,
00426                                  crSet.values_[j]))
00427       }
00428    }
00429 
00430    for(i=0; i<data_->numCRPenSets_; i++) {
00431       CRSet& crSet = data_->crPenSets_[i];
00432 
00433       for(int j=0; j<1; j++) {
00434          CHK_ERR(fei_->loadCRPen(crSet.crID_,
00435                                 crSet.numNodes_,
00436                                 crSet.nodeIDs_[j],
00437                                 crSet.fieldIDs_,
00438                                 crSet.weights_,
00439                                 crSet.values_[j],
00440                                 crSet.penValues_[j]))
00441       }
00442    }
00443   return(0);
00444 }
00445 
00446 //----------------------------------------------------------------------------
00447 int FEI_tester::aggregateLoadPhase()
00448 {
00449    int i;
00450 
00451    for(i=0; i<numMatrices; i++) {
00452       CHK_ERR(fei_->setCurrentMatrix(matrixIDs[i]))
00453 
00454       for(int j=0; j<data_->numElemBlocks_; j++) {
00455          ElemBlock& block = data_->elemBlocks_[j];
00456 
00457          for(int el=0; el<block.numElements_; el++) {
00458 
00459             CHK_ERR(fei_->sumInElemMatrix(block.blockID_,
00460                                    block.elemIDs_[el],
00461                                    block.elemConn_[el],
00462                                    block.elemStiff_[el],
00463                                    block.elemFormat_))
00464          }
00465       }
00466    }
00467 
00468    for(i=0; i<numRHSs; i++) {
00469       CHK_ERR(fei_->setCurrentRHS(rhsIDs[i]))
00470 
00471       for(int j=0; j<data_->numElemBlocks_; j++) {
00472          ElemBlock& block = data_->elemBlocks_[j];
00473 
00474          for(int el=0; el<block.numElements_; el++) {
00475             CHK_ERR(fei_->sumInElemRHS(block.blockID_,
00476                                    block.elemIDs_[el],
00477                                    block.elemConn_[el],
00478                                    block.elemLoad_[el]))
00479          }
00480       }
00481    }
00482 
00483    //*************** Boundary Condition Node Sets *************************
00484 
00485    for(i=0; i<data_->numBCNodeSets_; i++) {
00486       BCNodeSet& bcSet = data_->bcNodeSets_[i];
00487 
00488       CHK_ERR(fei_->loadNodeBCs(bcSet.numNodes_,
00489                      bcSet.nodeIDs_,
00490                      bcSet.fieldID_,
00491                      bcSet.offsetsIntoField_,
00492                      bcSet.prescribed_values_))
00493    }
00494 
00495    double* matScalars = new double[numMatrices];
00496    for(i=0; i<numMatrices; i++) {
00497       matScalars[i] = 1.0;
00498    }
00499 
00500    int rhsScaleID = rhsIDs[0];
00501    double rhsScalar = 1.0;
00502 
00503    CHK_ERR(fei_->setMatScalars(numMatrices, matrixIDs, matScalars))
00504    CHK_ERR(fei_->setRHSScalars(1, &rhsScaleID, &rhsScalar))
00505 
00506    delete [] matScalars;
00507   return(0);
00508 }
00509 
00510 //----------------------------------------------------------------------------
00511 int FEI_tester::exerciseResidualNorm()
00512 {
00513   std::string sname(data_->solverLibraryName_);
00514   if (sname == "TEST_LSC") {
00515     return( 0 );
00516   }
00517   //FEI_COUT << "numFields: " << data_->numFields_ << FEI_ENDL;
00518   double* norms = new double[data_->numFields_];
00519   int *fields = new int[data_->numFields_];
00520   for(int i=0; i<data_->numFields_; ++i) {
00521     fields[i] = data_->fieldIDs_[i];
00522   }
00523 
00524   CHK_ERR( fei_->residualNorm(1, data_->numFields_, fields, norms) );
00525   //int n;
00526   //for(n=0; n<data_->numFields_; n++) {
00527   //  FEI_COUT << " field["<<n<<"]: " << fields[n]
00528   //    << ", 1-norm: " << norms[n] << FEI_ENDL;
00529   //}
00530 
00531   delete [] fields;
00532   delete [] norms;
00533   return(0);
00534 }
00535 
00536 //----------------------------------------------------------------------------
00537 int FEI_tester::exercisePutFunctions()
00538 {
00539    //let's exercise the putNodalFieldData function.
00540 
00541    int numNodes;
00542    CHK_ERR( fei_->getNumLocalNodes(numNodes) );
00543    std::vector<int> nodeIDs(numNodes);
00544    int* nodeIDsPtr = &nodeIDs[0];
00545    int checkNumNodes;
00546    CHK_ERR( fei_->getLocalNodeIDList(checkNumNodes, nodeIDsPtr, numNodes) );
00547 
00548    for(int i=0; i<data_->numFields_; ++i) {
00549      int fieldID = data_->fieldIDs_[i];
00550      int fieldSize = data_->fieldSizes_[i];
00551      std::vector<double> data(numNodes*fieldSize, 0.0001);
00552 
00553      CHK_ERR( fei_->putNodalFieldData(fieldID, numNodes, nodeIDsPtr,
00554               &data[0]) );
00555    }
00556 
00557    return(0);
00558 }
00559 
00560 //----------------------------------------------------------------------------
00561 int FEI_tester::save_block_node_soln(DataReader& data, FEI& fei,
00562              const char* solnFileName, int numProcs,
00563              int localProc, int solveCounter)
00564 {
00565   (void)solveCounter;
00566    int i;
00567 
00568    int maxNumEqnsPerNode = 0;
00569    for(i=0; i<data.numFields_; ++i) {
00570      maxNumEqnsPerNode += data.fieldSizes_[i];
00571    }
00572 
00573    std::vector<double> soln(maxNumEqnsPerNode);
00574 
00575    int numNodes;
00576    CHK_ERR( fei.getNumLocalNodes(numNodes) );
00577 
00578    std::vector<GlobalID> nodes(numNodes);
00579    int* nodesPtr = &nodes[0];
00580 
00581    int checkNumNodes;
00582    CHK_ERR( fei.getLocalNodeIDList( checkNumNodes, nodesPtr, numNodes) );
00583 
00584    if (checkNumNodes != numNodes) {
00585      ERReturn(-1);
00586    }
00587 
00588    FEI_OSTRINGSTREAM fileName;
00589    fileName << solnFileName<<".node."<<solveCounter<<"."<<numProcs<<"."<<localProc;
00590    FEI_OFSTREAM outfile(fileName.str().c_str());
00591 
00592    if (!outfile || outfile.bad()) {
00593       FEI_COUT << "ERROR opening solution output file " << fileName.str() << FEI_ENDL;
00594       return(1);
00595    }
00596 
00597    outfile.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00598 
00599    std::vector<int> offsets(2);
00600 
00601    for(i=0; i<numNodes; ++i) {
00602      CHK_ERR( fei.getNodalSolution(1, &(nodesPtr[i]),
00603            &offsets[0], &soln[0]) );
00604 
00605      int numDOF = offsets[1];
00606 
00607      outfile << nodesPtr[i] << " " << numDOF << FEI_ENDL;
00608      for(int j=0; j<numDOF; j++) {
00609        outfile << soln[j] << " ";
00610      }
00611      outfile << FEI_ENDL;
00612    }
00613 
00614    return(0);
00615 }
00616 
00617 //----------------------------------------------------------------------------
00618 int FEI_tester::save_block_elem_soln(DataReader& data, FEI& fei,
00619              const char* solnFileName,
00620              int numProcs, int localProc,
00621              int solveCounter)
00622 {
00623    int returnValue = 0;
00624    FEI_OSTRINGSTREAM fileName;
00625    fileName << solnFileName<<".elem."<<solveCounter<<"."<<numProcs<<"."<<localProc;
00626    FEI_OFSTREAM outfile(fileName.str().c_str());
00627 
00628    if (!outfile || outfile.bad()) {
00629       FEI_COUT << "ERROR opening elem-solution output file " << fileName.str() << FEI_ENDL;
00630       return(1);
00631    }
00632 
00633    for(int i=0; i<data.numElemBlocks_; i++) {
00634       if (returnValue != 0) break;
00635 
00636       ElemBlock& eb = data.elemBlocks_[i];
00637 
00638       GlobalID blockID = eb.blockID_;
00639       int numElems;
00640       CHK_ERR( fei.getNumBlockElements(blockID, numElems))
00641 
00642       int dofPerElem;
00643       CHK_ERR( fei.getNumBlockElemDOF(blockID, dofPerElem))
00644       int totalNumElemDOF = numElems*dofPerElem;
00645 
00646       if (totalNumElemDOF < 1) {
00647   continue;
00648       }
00649 
00650       GlobalID* elemIDs = new GlobalID[numElems];
00651       if (elemIDs==NULL) return(-1);
00652 
00653       int err = fei.getBlockElemIDList(blockID, numElems, elemIDs);
00654       if (err) returnValue = 1;
00655 
00656       int* offsets = new int[numElems+1];
00657       if (offsets == NULL) return(-1);
00658 
00659       if (totalNumElemDOF > 0) {
00660          double* solnValues = new double[totalNumElemDOF];
00661          if (solnValues == NULL) return(-1);
00662 
00663          err = fei.getBlockElemSolution(blockID, numElems, elemIDs,
00664                                          dofPerElem, solnValues);
00665          if (err) returnValue = 1;
00666 
00667          if (!err) {
00668             for(int j=0; j<numElems; j++) {
00669 
00670                outfile << (int)elemIDs[j] << " " << dofPerElem << FEI_ENDL << "  ";
00671                for(int k=0; k<dofPerElem; k++) {
00672                   outfile << solnValues[j*dofPerElem + k] << " ";
00673                }
00674                outfile << FEI_ENDL;
00675             }
00676          }
00677 
00678          delete [] solnValues;
00679       }
00680 
00681       delete [] elemIDs;
00682       delete [] offsets;
00683    }
00684 
00685    return(0);
00686 }
00687 
00688 //----------------------------------------------------------------------------
00689 int FEI_tester::save_multiplier_soln(DataReader& data, FEI& fei,
00690        const char* solnFileName,
00691                           int numProcs, int localProc, int solveCounter)
00692 {
00693    int numCRs = 0;
00694 
00695    CHK_ERR( fei.getNumCRMultipliers(numCRs) );
00696 
00697    int* globalNumCRs = new int[numProcs];
00698 #ifndef FEI_SER
00699    if (MPI_Allgather(&numCRs, 1, MPI_INT, globalNumCRs, 1, MPI_INT,
00700      comm_) != MPI_SUCCESS) {
00701      ERReturn(-1);
00702    }
00703 #endif
00704 
00705    int localCRStart = 0;
00706    for(int p=0; p<localProc; p++) localCRStart += globalNumCRs[p];
00707 
00708    delete [] globalNumCRs;
00709 
00710    FEI_OSTRINGSTREAM fileName;
00711    fileName << solnFileName<<".mult."<<solveCounter<<"."<<numProcs<<"."<<localProc;
00712    FEI_OFSTREAM outfile(fileName.str().c_str());
00713 
00714    if (!outfile || outfile.bad()) {
00715       FEI_COUT << "ERROR opening mult-solution output file " << fileName.str() << FEI_ENDL;
00716       return(-1);
00717    }
00718 
00719    int* CRIDs = numCRs > 0 ? new int[numCRs] : NULL;
00720    double* results = numCRs > 0 ? new double[numCRs] : NULL;
00721 
00722    if (numCRs > 0 && (CRIDs==NULL || results==NULL)) {
00723      ERReturn(-1);
00724    }
00725 
00726    if (numCRs < 1) {
00727      return(0);
00728    }
00729 
00730    CHK_ERR( fei.getCRMultIDList(numCRs, CRIDs) );
00731 
00732    std::string sname(data_->solverLibraryName_);
00733    if (sname == "FETI") {
00734      for(int ii=0; ii<numCRs; ++ii) results[ii] = -999.99;
00735    }
00736    else {
00737      CHK_ERR( fei.getCRMultipliers(numCRs, CRIDs, results));
00738    }
00739 
00740    for(int i=0; i<numCRs; i++) {
00741       outfile << localCRStart++ << " " << 1 << FEI_ENDL;
00742 
00743       outfile << "   " << results[i] << FEI_ENDL;
00744    }
00745 
00746    delete [] CRIDs;
00747    delete [] results;
00748 
00749    return(0);
00750 }
00751 
00752 //----------------------------------------------------------------------------
00753 int FEI_tester::lsc_matrix_check()
00754 {
00755    if (localProc_ == 0) {
00756      char* current_dir = NULL;
00757      CHK_ERR( fei_test_utils::dirname(data_->solnFileName_.c_str(), current_dir));
00758 
00759      FEI_OSTRINGSTREAM solnMtxName;
00760      solnMtxName<< current_dir<<"/A_TLSC.mtx";
00761      fei::FillableMat solnMtx, checkMtx;
00762      CHK_ERR( SolnCheck::readMatrix(solnMtxName.str().c_str(), numProcs_, solnMtx) );
00763      CHK_ERR( SolnCheck::readMatrix(data_->checkFileName_.c_str(), numProcs_, checkMtx) );
00764      int err = SolnCheck::compareMatrices(solnMtx, checkMtx);
00765      delete [] current_dir;
00766      if (err == 0) {
00767        FEI_COUT << "Utst_fei_lsc: TEST PASSED" << FEI_ENDL;
00768      }
00769      else {
00770        FEI_COUT << "Utst_fei_lsc: TEST FAILED" << FEI_ENDL;
00771        return(-1);
00772      }
00773    }
00774 
00775    return(0);
00776 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends