snl_fei_tester.cpp

00001 /*--------------------------------------------------------------------*/
00002 /*    Copyright 2005 Sandia Corporation.                              */
00003 /*    Under the terms of Contract DE-AC04-94AL85000, there is a       */
00004 /*    non-exclusive license for use of this work by or on behalf      */
00005 /*    of the U.S. Government.  Export of this program may require     */
00006 /*    a license from the United States Government.                    */
00007 /*--------------------------------------------------------------------*/
00008 
00009 #include <fei_macros.hpp>
00010 
00011 #include <test_utils/snl_fei_tester.hpp>
00012 
00013 #include <fei_LinearSystemCore.hpp>
00014 #include <snl_fei_ArrayUtils.hpp>
00015 #include <test_utils/LibraryFactory.hpp>
00016 
00017 #include <fei_base.hpp>
00018 
00019 #ifdef HAVE_FEI_FETI
00020 #include <FETI_DP_FiniteElementData.h>
00021 #endif
00022 
00023 #include <test_utils/DataReader.hpp>
00024 #include <test_utils/SolnCheck.hpp>
00025 
00026 #undef fei_file
00027 #define fei_file "snl_fei_tester.cpp"
00028 #include <fei_ErrMacros.hpp>
00029 
00030 //----------------------------------------------------------------------------
00031 snl_fei_tester::snl_fei_tester(fei::SharedPtr<DataReader> data_reader,
00032                                MPI_Comm comm, int localProc, int numProcs)
00033   : comm_(comm),
00034     factory_(),
00035     vecSpace_(),
00036     matrixGraph_(),
00037     A_(),
00038     x_(),
00039     b_(),
00040     linSys_(NULL),
00041     linSysCore_(NULL),
00042     feData_(NULL),
00043     data_(data_reader),
00044     idTypes_(),
00045     numPatterns_(0),
00046     localProc_(localProc),
00047     numProcs_(numProcs)
00048 {
00049 }
00050 
00051 //----------------------------------------------------------------------------
00052 snl_fei_tester::~snl_fei_tester()
00053 {
00054   delete linSysCore_;
00055   delete feData_;
00056 }
00057 
00058 //----------------------------------------------------------------------------
00059 int snl_fei_tester::testInitialization()
00060 {
00061   if (factory_.get() == NULL) {
00062     try {
00063       factory_ = fei::create_fei_Factory(comm_, data_->solverLibraryName_.c_str());
00064     }
00065     catch (std::runtime_error& exc) {
00066       FEI_CERR << exc.what()<<FEI_ENDL;
00067       return(1);
00068     }
00069     if (factory_.get() == NULL) {
00070       ERReturn(-1);
00071     }
00072   }
00073 
00074   std::vector<std::string> stdstrings;
00075   fei::utils::char_ptrs_to_strings(data_->numParams_, data_->paramStrings_,
00076                                    stdstrings);
00077   fei::ParameterSet paramset;
00078   fei::utils::parse_strings(stdstrings, " ", paramset);
00079 
00080   if (!path_.empty()) {
00081     paramset.add(fei::Param("debugOutput", path_.c_str()));
00082   }
00083 
00084   factory_->parameters(paramset);
00085 
00086   vecSpace_ = factory_->createVectorSpace(comm_, NULL);
00087 
00088   vecSpace_->setParameters(paramset);
00089 
00090   defineFieldsAndIDTypes();
00091 
00092   fei::SharedPtr<fei::VectorSpace> dummy;
00093   matrixGraph_ = factory_->createMatrixGraph(vecSpace_, dummy, NULL);
00094 
00095   matrixGraph_->setParameters(paramset);
00096 
00097   CHK_ERR( initElemBlocks() );
00098 
00099   CHK_ERR( initConstraints() );
00100 
00101   int i;
00102   for(i=0; i<data_->numSharedNodeSets_; ++i) {
00103     CommNodeSet& nodeSet = data_->sharedNodeSets_[i];
00104 
00105     CHK_ERR( vecSpace_->initSharedIDs(nodeSet.numNodes_,
00106                                        idTypes_[nodeTypeOffset_],
00107                                        nodeSet.nodeIDs_,
00108                                        nodeSet.procsPerNode_,
00109                                        nodeSet.procs_) );
00110   }
00111 
00112   CHK_ERR( matrixGraph_->initComplete() );
00113 
00114   return(0);
00115 }
00116 
00117 //----------------------------------------------------------------------------
00118 void snl_fei_tester::dumpMatrixFiles()
00119 {
00120   FEI_OSTRINGSTREAM osstr;
00121   osstr << "A_" << A_->typeName() << ".np"<<numProcs_;
00122   std::string str = osstr.str();
00123   A_->writeToFile(str.c_str());
00124 }
00125 
00126 //----------------------------------------------------------------------------
00127 void snl_fei_tester::setParameter(const char* param)
00128 {
00129   std::vector<std::string> stdstrings;
00130   fei::utils::char_ptrs_to_strings(1, &param, stdstrings);
00131   fei::ParameterSet paramset;
00132   fei::utils::parse_strings(stdstrings, " ", paramset);
00133   factory_->parameters(paramset);
00134   vecSpace_->setParameters(paramset);
00135   matrixGraph_->setParameters(paramset);
00136 
00137   linSys_->parameters(1, &param);
00138   A_->parameters(paramset);
00139 }
00140 
00141 //----------------------------------------------------------------------------
00142 int snl_fei_tester::testLoading()
00143 {
00144   linSys_ = factory_->createLinearSystem(matrixGraph_);
00145 
00146   A_ = factory_->createMatrix(matrixGraph_);
00147   x_ = factory_->createVector(matrixGraph_, true);
00148   b_ = factory_->createVector(matrixGraph_);
00149 
00150   matrixGraph_->setIndicesMode(fei::MatrixGraph::POINT_ENTRY_GRAPH);
00151 
00152   CHK_ERR( linSys_->parameters(data_->numParams_, data_->paramStrings_) );
00153 
00154   std::vector<std::string> stdstrings;
00155   fei::utils::char_ptrs_to_strings(data_->numParams_, data_->paramStrings_, stdstrings);
00156   fei::ParameterSet paramset;
00157   fei::utils::parse_strings(stdstrings, " ", paramset);
00158   CHK_ERR( A_->parameters(paramset) );
00159 
00160   linSys_->setMatrix(A_);
00161   linSys_->setRHS(b_);
00162   linSys_->setSolutionVector(x_);
00163 
00164   CHK_ERR( A_->putScalar(0.0) );
00165   CHK_ERR( b_->putScalar(0.0) );
00166 
00167   matrixGraph_->createSlaveMatrices();
00168 
00169   CHK_ERR( loadElemBlocks() );
00170 
00171   CHK_ERR( loadConstraints() );
00172 
00173   int i;
00174   for(i=0; i<data_->numBCNodeSets_; ++i) {
00175     BCNodeSet& bcSet = data_->bcNodeSets_[i];
00176     int fieldSize = data_->getFieldSize(bcSet.fieldID_);
00177     if (fieldSize < 1) {
00178       continue;
00179     }
00180 
00181     CHK_ERR( linSys_->loadEssentialBCs(bcSet.numNodes_,
00182                                        bcSet.nodeIDs_,
00183                                        idTypes_[nodeTypeOffset_],
00184                                        bcSet.fieldID_,
00185                                        bcSet.offsetsIntoField_,
00186                                        bcSet.prescribed_values_) );
00187   }
00188 
00189   CHK_ERR( linSys_->loadComplete() );
00190 
00191   return(0);
00192 }
00193 
00194 //----------------------------------------------------------------------------
00195 int snl_fei_tester::testSolve()
00196 {
00197   fei::SharedPtr<fei::Solver> solver = factory_->createSolver();
00198 
00199   std::vector<std::string> stdstrings;
00200   fei::utils::char_ptrs_to_strings(data_->numParams_, data_->paramStrings_, stdstrings);
00201   fei::ParameterSet paramset;
00202   fei::utils::parse_strings(stdstrings," ",paramset);
00203 
00204   int status, itersTaken = 0;
00205   CHK_ERR( solver->solve(linSys_.get(),
00206                          NULL, //preconditioningMatrix
00207                          paramset, itersTaken, status) );
00208 
00209   CHK_ERR( x_->scatterToOverlap() );
00210 
00211   return(0);
00212 }
00213 
00214 //----------------------------------------------------------------------------
00215 int snl_fei_tester::testCheckResult()
00216 {
00217   CHK_ERR( save_block_node_soln(*data_, x_.get(), data_->solnFileName_.c_str(),
00218                                 numProcs_, localProc_, 1));
00219 
00220   CHK_ERR( save_block_elem_soln(*data_, x_.get(), data_->solnFileName_.c_str(),
00221                                 numProcs_, localProc_, 1));
00222 
00223   CHK_ERR( save_multiplier_soln(*data_, x_.get(), data_->solnFileName_.c_str(),
00224                                 numProcs_, localProc_, 1));
00225 
00226   int err = SolnCheck::checkSolution(localProc_, numProcs_, data_->solnFileName_.c_str(),
00227                                      data_->checkFileName_.c_str(), "node", 1);
00228 
00229   err += SolnCheck::checkSolution(localProc_, numProcs_, data_->solnFileName_.c_str(),
00230                                   data_->checkFileName_.c_str(), "elem", 1);
00231 
00232   err += SolnCheck::checkSolution(localProc_, numProcs_, data_->solnFileName_.c_str(),
00233                                   data_->checkFileName_.c_str(), "mult", 1);
00234   int globalErr = err;
00235 #ifndef FEI_SER
00236   if (MPI_SUCCESS != MPI_Allreduce(&err, &globalErr, 1, MPI_INT, MPI_SUM,
00237                                    comm_)) return(-1);
00238 #endif
00239   if (globalErr != 0) return(-1);
00240   return(0);
00241 }
00242 
00243 //----------------------------------------------------------------------------
00244 void snl_fei_tester::defineFieldsAndIDTypes()
00245 {
00246   vecSpace_->defineFields(data_->numFields_, data_->fieldIDs_, data_->fieldSizes_);
00247 
00248   //nodeIDType == 0
00249   idTypes_.push_back(0);
00250 
00251   //constraintIDType == 1
00252   idTypes_.push_back(1);
00253 
00254   //elemDofIDType == 2
00255   idTypes_.push_back(2);
00256 
00257   vecSpace_->defineIDTypes(idTypes_.size(), &idTypes_[0] );
00258 
00259   nodeTypeOffset_ = 0;
00260   constraintTypeOffset_ = 1;
00261   elemTypeOffset_ = 2;
00262 }
00263 
00264 //----------------------------------------------------------------------------
00265 int snl_fei_tester::initElemBlocks()
00266 {
00267   for(int i=0; i<data_->numElemBlocks_; ++i) {
00268     ElemBlock& eb = data_->elemBlocks_[i];
00269 
00270     int patternID;
00271     definePattern(eb, patternID);
00272 
00273     CHK_ERR( matrixGraph_->initConnectivityBlock(eb.blockID_,
00274                                                eb.numElements_,
00275                                                patternID) );
00276 
00277     for(int j=0; j<eb.numElements_; ++j) {
00278       std::vector<int> conn(eb.numNodesPerElement_);
00279       for(int ii=0; ii<eb.numNodesPerElement_; ++ii) {
00280         conn[ii] = eb.elemConn_[j][ii];
00281       }
00282 
00283       CHK_ERR( matrixGraph_->initConnectivity(eb.blockID_,
00284                                               eb.elemIDs_[j],
00285                                               &conn[0]) );
00286     }
00287   }
00288 
00289   return(0);
00290 }
00291 
00292 //----------------------------------------------------------------------------
00293 int snl_fei_tester::loadElemBlocks()
00294 {
00295   int i;
00296   for(i=0; i<data_->numElemBlocks_; ++i) {
00297     ElemBlock& eb = data_->elemBlocks_[i];
00298 
00299     if (eb.numElements_ < 1) {
00300       continue;
00301     }
00302 
00303     int numIndices = matrixGraph_->getConnectivityNumIndices(eb.blockID_);
00304 
00305     std::vector<int> indices(numIndices);
00306 
00307     for(int j=0; j<eb.numElements_; ++j) {
00308       int checkNum;
00309       CHK_ERR( matrixGraph_->getConnectivityIndices(eb.blockID_,
00310                                                     eb.elemIDs_[j],
00311                                                     numIndices,
00312                                                     &indices[0],
00313                                                     checkNum) );
00314       if (numIndices != checkNum) {
00315         ERReturn(-1);
00316       }
00317 
00318       CHK_ERR( A_->sumIn(eb.blockID_, eb.elemIDs_[j],
00319                          eb.elemStiff_[j]) );
00320 
00321       CHK_ERR( b_->sumIn(numIndices, &indices[0],
00322                          eb.elemLoad_[j], 0) );
00323     }
00324   }
00325 
00326   return(0);
00327 }
00328 
00329 //----------------------------------------------------------------------------
00330 int snl_fei_tester::initConstraints()
00331 {
00332   std::vector<int> idTypes;
00333   int constraintID = localProc_*100000;
00334   int i;
00335   for(i=0; i<data_->numCRMultSets_; ++i) {
00336     CRSet& crSet = data_->crMultSets_[i];
00337 
00338     for(int j=0; j<1; ++j) {
00339       idTypes.assign(crSet.numNodes_, idTypes_[nodeTypeOffset_]);
00340 
00341       crSet.crID_ = constraintID++;
00342       int constraintIDType = idTypes_[constraintTypeOffset_];
00343       CHK_ERR( matrixGraph_->initLagrangeConstraint(crSet.crID_,
00344                                                   constraintIDType,
00345                                                   crSet.numNodes_,
00346                                                   &idTypes[0],
00347                                                   crSet.nodeIDs_[j],
00348                                                   crSet.fieldIDs_) );
00349     }
00350   }
00351 
00352   for(i=0; i<data_->numCRPenSets_; ++i) {
00353     CRSet& crSet = data_->crPenSets_[i];
00354 
00355     for(int j=0; j<1; ++j) {
00356       idTypes.assign(crSet.numNodes_, idTypes_[nodeTypeOffset_]);
00357 
00358       crSet.crID_ = constraintID++;
00359       int constraintIDType = idTypes_[constraintTypeOffset_];
00360       CHK_ERR( matrixGraph_->initPenaltyConstraint(crSet.crID_,
00361                                                   constraintIDType,
00362                                                   crSet.numNodes_,
00363                                                   &idTypes[0],
00364                                                   crSet.nodeIDs_[j],
00365                                                   crSet.fieldIDs_) );
00366     }
00367   }
00368 
00369   std::map<int,int> fieldDB;
00370   for(i=0; i<data_->numFields_; ++i) {
00371     fieldDB.insert(std::pair<int,int>(data_->fieldIDs_[i], data_->fieldSizes_[i]));
00372   }
00373 
00374   std::vector<int> nodeIDs;
00375   std::vector<int> fieldIDs;
00376   std::vector<double> weights;
00377 
00378   for(i=0; i<data_->numSlaveVars_; i++) {
00379     int ii;
00380     CRSet& crSet = data_->slaveVars_[i];
00381 
00382     nodeIDs.resize(crSet.numNodes_+1);
00383     nodeIDs[0] = crSet.slaveNodeID_;
00384     fieldIDs.resize(0);
00385     fieldIDs.push_back(crSet.slaveFieldID_);
00386 
00387     for(ii=0; ii<crSet.numNodes_; ++ii) {
00388       nodeIDs[ii+1] = crSet.nodeIDs_[0][ii];
00389       fieldIDs.push_back(crSet.fieldIDs_[ii]);
00390     }
00391 
00392     idTypes.assign(crSet.numNodes_+1, idTypes_[nodeTypeOffset_]);
00393 
00394     int fieldSize = fieldDB[crSet.slaveFieldID_];
00395     weights.resize(0);
00396     for(ii=0; ii<fieldSize; ++ii) weights.push_back(0.0);
00397     weights[crSet.slaveOffset_] = -1.0;
00398     int offset = 0;
00399     for(ii=0; ii<crSet.numNodes_; ++ii) {
00400       fieldSize = fieldDB[crSet.fieldIDs_[ii]];
00401       for(int jj=0; jj<fieldSize; ++jj) {
00402         weights.push_back(crSet.weights_[offset++]);
00403       }
00404     }
00405 
00406     CHK_ERR( matrixGraph_->initSlaveConstraint(crSet.numNodes_+1,
00407                                                &idTypes[0],
00408                                                &nodeIDs[0],
00409                                                &fieldIDs[0],
00410                                                0,
00411                                                crSet.slaveOffset_,
00412                                                &weights[0],
00413                                                crSet.values_[0]));
00414   }
00415 
00416   return(0);
00417 }
00418 
00419 //----------------------------------------------------------------------------
00420 int snl_fei_tester::loadConstraints()
00421 {
00422   int i;
00423   for(i=0; i<data_->numCRMultSets_; ++i) {
00424     CRSet& crSet = data_->crMultSets_[i];
00425 
00426     for(int j=0; j<1; ++j) {
00427       CHK_ERR( linSys_->loadLagrangeConstraint(crSet.crID_,
00428                                                crSet.weights_,
00429                                                crSet.values_[j]) );
00430     }
00431   }
00432 
00433   for(i=0; i<data_->numCRPenSets_; ++i) {
00434     CRSet& crSet = data_->crPenSets_[i];
00435 
00436     for(int j=0; j<1; ++j) {
00437       CHK_ERR( linSys_->loadPenaltyConstraint(crSet.crID_,
00438                                               crSet.weights_,
00439                                               crSet.penValues_[j],
00440                                               crSet.values_[j]) );
00441     }
00442   }
00443 
00444   return(0);
00445 }
00446 
00447 //----------------------------------------------------------------------------
00448 void snl_fei_tester::definePattern(ElemBlock& eb, int& patternID)
00449 {
00450   int i, j, numIDTypes = 1;
00451   numIDTypes += eb.numElemDOF_>0 ? 1 : 0;
00452 
00453   //find out how many nodal fields there are, total.
00454   std::vector<int> nodalFieldIDs;
00455   std::vector<int> flatFieldIDsArray;
00456   for(i=0; i<eb.numNodesPerElement_; ++i) {
00457     for(j=0; j<eb.numFieldsPerNode_[i]; ++j) {
00458       snl_fei::sortedListInsert(eb.nodalFieldIDs_[i][j], nodalFieldIDs);
00459       flatFieldIDsArray.push_back(eb.nodalFieldIDs_[i][j]);
00460     }
00461   }
00462 
00463   patternID = numPatterns_++;
00464 
00465   if (numIDTypes == 1 && nodalFieldIDs.size() == 1) {
00466     //This is a very simple pattern
00467     patternID = matrixGraph_->definePattern(eb.numNodesPerElement_,
00468                                 idTypes_[nodeTypeOffset_],
00469                                 nodalFieldIDs[0]);
00470   }
00471   else if (numIDTypes == 1) {
00472     std::vector<int> numFieldsPerID(eb.numNodesPerElement_);
00473 
00474     patternID = matrixGraph_->definePattern(eb.numNodesPerElement_,
00475                                 idTypes_[nodeTypeOffset_],
00476                                 eb.numFieldsPerNode_,
00477                                 &flatFieldIDsArray[0]);
00478   }
00479   else {
00480     std::vector<int> idTypes(eb.numNodesPerElement_+1, idTypes_[nodeTypeOffset_]);
00481     idTypes[idTypes.size()-1] = idTypes_[elemTypeOffset_];
00482     std::vector<int> numFieldsPerID(idTypes.size());
00483     std::vector<int> fieldIDs;
00484     for(i=0; i<eb.numNodesPerElement_; ++i) {
00485       numFieldsPerID[i] = eb.numFieldsPerNode_[i];
00486       for(j=0; j<eb.numFieldsPerNode_[i]; ++j) {
00487         fieldIDs.push_back(eb.nodalFieldIDs_[i][j]);
00488       }
00489     }
00490     numFieldsPerID[idTypes.size()-1] = eb.numElemDOF_;
00491     for(i=0; i<eb.numElemDOF_; ++i) {
00492       fieldIDs.push_back(eb.elemDOFFieldIDs_[i]);
00493     }
00494 
00495     patternID = matrixGraph_->definePattern(idTypes.size(),
00496                                 &idTypes[0],
00497                                 &numFieldsPerID[0],
00498                                 &fieldIDs[0]);
00499   }
00500 }
00501 
00502 //----------------------------------------------------------------------------
00503 int snl_fei_tester::save_block_node_soln(DataReader& data, fei::Vector* vec,
00504                                      const char* solnFileName, int numProcs,
00505                                      int localProc, int solveCounter)
00506 {
00507   (void)solveCounter;
00508 
00509   int numLocalNodes = vecSpace_->getNumOwnedAndSharedIDs(idTypes_[nodeTypeOffset_]);
00510 
00511   int* nodeList = new int[numLocalNodes];
00512 
00513   int checkNum = 0;
00514   int err = vecSpace_->getOwnedAndSharedIDs(idTypes_[nodeTypeOffset_],
00515                                     numLocalNodes, nodeList, checkNum);
00516   if (err != 0) {
00517     ERReturn(-1);
00518   }
00519 
00520   FEI_OSTRINGSTREAM fileName;
00521   fileName<< solnFileName<<".node."<<solveCounter<<"."<<numProcs<<"."<<localProc;
00522   FEI_OFSTREAM outfile(fileName.str().c_str());
00523 
00524   if (!outfile || outfile.bad()) {
00525     FEI_CERR << "ERROR opening solution output file " << fileName << FEI_ENDL;
00526     return(-1);
00527   }
00528 
00529   outfile.setf(IOS_SCIENTIFIC, IOS_FLOATFIELD);
00530 
00531   std::vector<double> solnData;
00532   std::vector<int> fieldList;
00533 
00534   for(int i=0; i<numLocalNodes; i++) {
00535     int idType = idTypes_[nodeTypeOffset_];
00536     int ID = nodeList[i];
00537 
00538     int numDOF = vecSpace_->getNumDegreesOfFreedom(idType, ID);
00539     int numFields = vecSpace_->getNumFields(idType, ID);
00540     solnData.resize(numDOF);
00541     vecSpace_->getFields(idType, ID, fieldList);
00542 
00543     outfile << ID << " " << numDOF << FEI_ENDL;
00544     for(int j=0; j<numFields; ++j) {
00545       int fieldSize = vecSpace_->getFieldSize(fieldList[j]);
00546 
00547       CHK_ERR( vec->copyOutFieldData(fieldList[j], idType,
00548                                      1, &ID, &solnData[0]) );
00549 
00550       for(int k=0; k<fieldSize; ++k) {
00551         outfile << solnData[k] << " ";
00552       }
00553     }
00554     outfile << FEI_ENDL;
00555   }
00556 
00557   delete [] nodeList;
00558 
00559   return(0);
00560 }
00561 
00562 //----------------------------------------------------------------------------
00563 int snl_fei_tester::save_block_elem_soln(DataReader& data, fei::Vector* vec,
00564                                          const char* solnFileName,
00565                                          int numProcs,
00566                                          int localProc, int solveCounter)
00567 {
00568   (void)solveCounter;
00569 
00570   int numLocalElems = vecSpace_->getNumOwnedAndSharedIDs(idTypes_[elemTypeOffset_]);
00571 
00572   int* elemList = new int[numLocalElems];
00573 
00574   int checkNum = 0;
00575   int err = vecSpace_->getOwnedAndSharedIDs(idTypes_[elemTypeOffset_],
00576                                     numLocalElems, elemList, checkNum);
00577   if (err != 0) {
00578     ERReturn(-1);
00579   }
00580 
00581   FEI_OSTRINGSTREAM fileName;
00582   fileName<< solnFileName<<".elem."<<solveCounter<<"."<<numProcs<<"."<<localProc;
00583   FEI_OFSTREAM outfile(fileName.str().c_str());
00584 
00585   if (!outfile || outfile.bad()) {
00586     FEI_CERR << "ERROR opening solution output file " << fileName << FEI_ENDL;
00587     return(-1);
00588   }
00589 
00590   std::vector<double> solnData;
00591   std::vector<int> fieldList;
00592 
00593   for(int i=0; i<numLocalElems; i++) {
00594     int idType = idTypes_[elemTypeOffset_];
00595     int ID = elemList[i];
00596 
00597     int numDOF = vecSpace_->getNumDegreesOfFreedom(idType, ID);
00598     int numFields = vecSpace_->getNumFields(idType, ID);
00599     solnData.resize(numDOF);
00600     vecSpace_->getFields(idType, ID, fieldList);
00601 
00602     outfile << ID << " " << numDOF << FEI_ENDL;
00603     for(int j=0; j<numFields; ++j) {
00604       int fieldSize = vecSpace_->getFieldSize(fieldList[j]);
00605 
00606       CHK_ERR( vec->copyOutFieldData(fieldList[j], idType,
00607                                      1, &ID, &solnData[0]) );
00608 
00609       for(int k=0; k<fieldSize; ++k) {
00610         outfile << solnData[k] << " ";
00611       }
00612     }
00613     outfile << FEI_ENDL;
00614   }
00615 
00616   delete [] elemList;
00617 
00618   return(0);
00619 }
00620 
00621 //----------------------------------------------------------------------------
00622 int snl_fei_tester::save_multiplier_soln(DataReader& data, fei::Vector* vec,
00623                                          const char* solnFileName,
00624                                          int numProcs, int localProc,
00625                                          int solveCounter)
00626 {
00627   (void)solveCounter;
00628 
00629   int numLocalCRs = vecSpace_->getNumOwnedAndSharedIDs(idTypes_[constraintTypeOffset_]);
00630 
00631    int* globalNumCRs = new int[numProcs];
00632 #ifndef FEI_SER
00633    if (MPI_Allgather(&numLocalCRs, 1, MPI_INT, globalNumCRs, 1, MPI_INT,
00634                  comm_) != MPI_SUCCESS) {
00635      ERReturn(-1);
00636    }
00637 #endif
00638 
00639    int localCRStart = 0;
00640 #ifndef FEI_SER
00641    for(int p=0; p<localProc; p++) localCRStart += globalNumCRs[p];
00642 #endif
00643 
00644    delete [] globalNumCRs;
00645 
00646   std::vector<int> crList(numLocalCRs);
00647 
00648   int checkNum = 0;
00649   int err = vecSpace_->getOwnedAndSharedIDs(
00650     idTypes_[constraintTypeOffset_], numLocalCRs,
00651     numLocalCRs ? &crList[0] : 0, checkNum);
00652   if (err != 0) {
00653     ERReturn(-1);
00654   }
00655 
00656   FEI_OSTRINGSTREAM fileName;
00657   fileName<< solnFileName<<".mult."<<solveCounter<<"."<<numProcs<<"."<<localProc;
00658   FEI_OFSTREAM outfile(fileName.str().c_str());
00659 
00660   if (!outfile || outfile.bad()) {
00661     FEI_CERR << "ERROR opening solution output file " << fileName << FEI_ENDL;
00662     return(-1);
00663   }
00664 
00665   std::vector<double> solnData;
00666   std::vector<int> fieldList;
00667 
00668   for(int i=0; i<numLocalCRs; i++) {
00669     int idType = idTypes_[constraintTypeOffset_];
00670     int ID = crList[i];
00671 
00672     solnData.resize(1);
00673 
00674     outfile << localCRStart++ << " " << 1 << FEI_ENDL;
00675     for(int j=0; j<1; ++j) {
00676       int globalIndex = -1;
00677       CHK_ERR( vecSpace_->getGlobalIndex(idType, ID, globalIndex) );
00678 
00679       CHK_ERR( vec->copyOut(1, &globalIndex, &solnData[0]) );
00680 
00681       for(int k=0; k<1; ++k) {
00682         outfile << solnData[k] << " ";
00683       }
00684     }
00685     outfile << FEI_ENDL;
00686   }
00687 
00688   return(0);
00689 }

Generated on Wed May 12 21:30:42 2010 for FEI by  doxygen 1.4.7