FEI Version of the Day
HexBeam.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_macros.hpp>
00045 
00046 #include <test_utils/HexBeam.hpp>
00047 
00048 #undef fei_file
00049 #define fei_file "HexBeam.cpp"
00050 
00051 HexBeam::HexBeam(int W, int D, int DofPerNode,
00052      int decomp, int numProcs, int localProc)
00053   : W_(W),
00054     D_(D),
00055     decomp_(decomp),
00056     numProcs_(numProcs),
00057     localProc_(localProc),
00058     firstLocalElem_(0),
00059     firstLocalNode_(0),
00060     inErrorState_(false),
00061     nodesPerElem_(8),
00062     dofPerNode_(DofPerNode)
00063 {
00064   totalNumElems_ = W*W*D;
00065   totalNumNodes_ = (W+1)*(W+1)*(D+1);
00066   numElemsPerSlice_ = W*W;
00067   numNodesPerSlice_ = (W+1)*(W+1);
00068 
00069   numGlobalDOF_ = totalNumNodes_*dofPerNode_;
00070 
00071   numLocalSlices_ = D/numProcs;
00072   int remainder = D%numProcs;
00073 
00074   switch(decomp) {
00075   case HexBeam::OneD:
00076     if (D < numProcs) {
00077       throw std::runtime_error("HexBeam: size D must be greater or equal num-procs.");
00078     }
00079     if (localProc < remainder) {
00080       ++numLocalSlices_;
00081     }
00082 
00083     localNumNodes_ = numNodesPerSlice_*(numLocalSlices_+1);
00084     localNumElems_ = numElemsPerSlice_*numLocalSlices_;
00085     numLocalDOF_ = localNumNodes_*dofPerNode_;
00086 
00087     if (localProc > 0) {
00088       firstLocalElem_ = localProc*numLocalSlices_*numElemsPerSlice_;
00089       firstLocalNode_ = localProc*numLocalSlices_*numNodesPerSlice_;
00090       if (remainder <= localProc && remainder > 0) {
00091   firstLocalElem_ += remainder*numElemsPerSlice_;
00092   firstLocalNode_ += remainder*numNodesPerSlice_;
00093       }
00094     }
00095 
00096     break;
00097 
00098   case HexBeam::TwoD:
00099   case HexBeam::ThreeD:
00100   default:
00101     fei::console_out() << "HexBeam: invalid decomp option: " << decomp
00102    <<" aborting." << FEI_ENDL;
00103     std::abort();
00104   }
00105 }
00106 
00107 HexBeam::~HexBeam()
00108 {
00109 }
00110 
00111 int HexBeam::getElemConnectivity(int elemID, int* nodeIDs)
00112 {
00113   if (elemID < firstLocalElem_ || elemID > firstLocalElem_+localNumElems_) {
00114     return(-1);
00115   }
00116 
00117   int whichGlobalSlice = elemID/numElemsPerSlice_;
00118   int elemX = elemID%W_;
00119   int elemY = (elemID%(W_*W_))/W_;
00120 
00121   int firstElemNode = whichGlobalSlice*numNodesPerSlice_
00122                      + elemY*(W_+1) + elemX;
00123 
00124   nodeIDs[0] = firstElemNode;
00125   nodeIDs[1] = firstElemNode+1;
00126   nodeIDs[2] = firstElemNode+W_+1;
00127   nodeIDs[3] = nodeIDs[2]+1;
00128 
00129   nodeIDs[4] = nodeIDs[0]+numNodesPerSlice_;
00130   nodeIDs[5] = nodeIDs[1]+numNodesPerSlice_;
00131   nodeIDs[6] = nodeIDs[2]+numNodesPerSlice_;
00132   nodeIDs[7] = nodeIDs[3]+numNodesPerSlice_;
00133 
00134   return(0);
00135 }
00136 
00137 int HexBeam::getElemStiffnessMatrix(int elemID, double* elemMat)
00138 {
00139   if (elemID < firstLocalElem_ || elemID > firstLocalElem_+localNumElems_) {
00140     return(-1);
00141   }
00142 
00143   int i, len = nodesPerElem_*dofPerNode_*nodesPerElem_*dofPerNode_;
00144 
00145   for(i=0; i<len; ++i) {
00146     elemMat[i] = 0.0;
00147   }
00148 
00149   //Should set up some semi-realistic stiffness-matrix coefficients here...
00150   //For now just use arbitrary numbers and set it up so the matrix won't be
00151   //too ill-conditioned. (This is intended for an assembly test more than
00152   //a solver test.)
00153 
00154   //Now set the diagonal to 4.0
00155   len = nodesPerElem_*dofPerNode_;
00156   for(i=0; i<len; ++i) {
00157     int offset = i*len+i;
00158     elemMat[offset] = 4.0;
00159   }
00160 
00161   //Now set some off-diagonals
00162   for(i=0; i<len; ++i) {
00163     int offset = i*len+i;
00164     if (i>1) {
00165       elemMat[offset-2] = -0.5;
00166     }
00167 
00168     if (i<len-2) {
00169       elemMat[offset+2] = -0.5;
00170     }
00171 
00172     if (i>3) {
00173       elemMat[offset-4] = -0.1;
00174     }
00175     if (i<len-4) {
00176       elemMat[offset+4] = -0.1;
00177     }
00178   }
00179 
00180   return(0);
00181 }
00182 
00183 int HexBeam::getElemLoadVector(int elemID, double* elemVec)
00184 {
00185   if (elemID < firstLocalElem_ || elemID > firstLocalElem_+localNumElems_) {
00186     return(-1);
00187   }
00188 
00189   int i, len = nodesPerElem_*dofPerNode_;
00190   for(i=0; i<len; ++i) {
00191     elemVec[i] = 1.0;
00192   }
00193 
00194   return(0);
00195 }
00196 
00197 int HexBeam::getNumBCNodes()
00198 {
00199   int numBCNodes = (numLocalSlices_+1)*(W_+1);
00200   return( numBCNodes );
00201 }
00202 
00203 int HexBeam::getBCNodes(int numNodes, int* nodeIDs)
00204 {
00205   if (numNodes != getNumBCNodes()) {
00206     return(-1);
00207   }
00208 
00209   int firstBCNode = firstLocalNode_ + W_;
00210 
00211   for(int i=0; i<numNodes; ++i) {
00212     nodeIDs[i] = firstBCNode + W_+1;
00213   }
00214 
00215   return(0);
00216 }
00217 
00218 int HexBeam::getBCValues(int numBCNodes, int* offsetsIntoField, double* vals)
00219 {
00220   if (numBCNodes != getNumBCNodes()) {
00221     return(-1);
00222   }
00223 
00224   for(int i=0; i<numBCNodes; ++i) {
00225     offsetsIntoField[i] = 0;
00226     vals[i] = 2.0;
00227   }
00228 
00229   return(0);
00230 }
00231 
00232 int HexBeam::getNumSharedNodes()
00233 {
00234   if (numProcs_ < 2) return(0);
00235 
00236   int numSharedNodes = numNodesPerSlice_;
00237   if (localProc_ > 0 && localProc_ < numProcs_-1) {
00238     numSharedNodes += numNodesPerSlice_;
00239   }
00240 
00241   return(numSharedNodes);
00242 }
00243 
00244 int HexBeam::getSharedNodes(int numSharedNodes,
00245           int*& sharedNodes,
00246           int*& numSharingProcsPerNode,
00247           int**& sharingProcs)
00248 {
00249   if (numProcs_ < 2) return(0);
00250 
00251   if (numSharedNodes != getNumSharedNodes()) {
00252     return(-1);
00253   }
00254 
00255   sharedNodes = new int[numSharedNodes];
00256   numSharingProcsPerNode = new int[numSharedNodes];
00257   sharingProcs = new int*[numSharedNodes];
00258   int* sharingProcVals = new int[numSharedNodes];
00259   if (sharedNodes == NULL || numSharingProcsPerNode == NULL ||
00260       sharingProcs == NULL || sharingProcVals == NULL) {
00261     return(-1);
00262   }
00263 
00264   int i;
00265   for(i=0; i<numSharedNodes; ++i) {
00266     numSharingProcsPerNode[i] = 1;
00267     sharingProcs[i] = &(sharingProcVals[i]);
00268   }
00269 
00270   int firstSharedNode = firstLocalNode_+numNodesPerSlice_*numLocalSlices_;
00271   int offset = 0;
00272 
00273   if (localProc_ < numProcs_ - 1) {
00274     for(i=0; i<numNodesPerSlice_; ++i) {
00275       sharedNodes[offset] = firstSharedNode+i;
00276       sharingProcs[offset++][0] = localProc_+1;
00277     }
00278   }
00279 
00280   firstSharedNode = firstLocalNode_;
00281   if (localProc_ > 0) {
00282     for(i=0; i<numNodesPerSlice_; ++i) {
00283       sharedNodes[offset] = firstSharedNode+i;
00284       sharingProcs[offset++][0] = localProc_-1;
00285     }
00286   }
00287 
00288   return(0);
00289 }
00290 
00291 namespace HexBeam_Functions {
00292 
00293 int print_cube_data(HexBeam& hexcube, int numProcs, int localProc)
00294 {
00295   FEI_COUT << localProc << ": num elems: " << hexcube.numLocalElems() << FEI_ENDL;
00296   int i;
00297   int* nodeIDs = new int[hexcube.numNodesPerElem()];
00298   int firstLocalElem = hexcube.firstLocalElem();
00299 
00300   for(i=0; i<hexcube.numLocalElems(); ++i) {
00301     hexcube.getElemConnectivity(firstLocalElem+i, nodeIDs);
00302     FEI_COUT << localProc << ": elem " << firstLocalElem+i << ", nodes: ";
00303     for(int j=0; j<hexcube.numNodesPerElem(); ++j) {
00304       FEI_COUT << nodeIDs[j] << " ";
00305     }
00306     FEI_COUT << FEI_ENDL;
00307   }
00308 
00309   delete [] nodeIDs;
00310 
00311   return(0);
00312 }
00313 
00314 
00315 int init_elem_connectivities(FEI* fei, HexBeam& hexcube)
00316 {
00317   int numLocalElems = hexcube.numLocalElems();
00318   int firstLocalElem = hexcube.firstLocalElem();
00319   int nodesPerElem = hexcube.numNodesPerElem();
00320   int fieldID = 0;
00321 
00322   int** fieldIDsTable = new int*[nodesPerElem];
00323   int* numFieldsPerNode = new int[nodesPerElem];
00324 
00325   for(int j=0; j<nodesPerElem; ++j) {
00326     numFieldsPerNode[j] = 1;
00327     fieldIDsTable[j] = new int[numFieldsPerNode[j]];
00328     for(int k=0; k<numFieldsPerNode[j]; ++k) {
00329       fieldIDsTable[j][k] = fieldID;
00330     }
00331   }
00332 
00333   int blockID = 0;
00334   CHK_ERR( fei->initElemBlock(blockID,
00335             numLocalElems,
00336             nodesPerElem,
00337             numFieldsPerNode,
00338             fieldIDsTable,
00339             0, // no element-centered degrees-of-freedom
00340             NULL, //null list of elem-dof fieldIDs
00341             FEI_NODE_MAJOR) );
00342 
00343 
00344   int* nodeIDs = new int[nodesPerElem];
00345   if (nodeIDs == NULL) return(-1);
00346 
00347   for(int i=0; i<numLocalElems; ++i) {
00348     CHK_ERR( hexcube.getElemConnectivity(firstLocalElem+i, nodeIDs) );
00349 
00350     CHK_ERR( fei->initElem(blockID, firstLocalElem+i, nodeIDs) );
00351   }
00352 
00353   delete [] nodeIDs;
00354   delete [] numFieldsPerNode;
00355   for(int jj=0; jj<nodesPerElem; ++jj) {
00356     delete [] fieldIDsTable[jj];
00357   }
00358   delete [] fieldIDsTable;
00359 
00360   return(0);
00361 }
00362 
00363 int init_shared_nodes(FEI* fei, HexBeam& hexcube)
00364 {
00365   int numSharedNodes = hexcube.getNumSharedNodes();
00366   if (numSharedNodes == 0) {
00367     return(0);
00368   }
00369 
00370   int* sharedNodes = NULL;
00371   int* numSharingProcsPerNode = NULL;
00372   int** sharingProcs = NULL;
00373   if (numSharedNodes > 0) {
00374     CHK_ERR( hexcube.getSharedNodes(numSharedNodes,
00375             sharedNodes, numSharingProcsPerNode,
00376             sharingProcs) );
00377   }
00378 
00379   CHK_ERR( fei->initSharedNodes(numSharedNodes, sharedNodes,
00380             numSharingProcsPerNode, sharingProcs) );
00381 
00382   delete [] sharedNodes;
00383   delete [] numSharingProcsPerNode;
00384   delete [] sharingProcs[0];
00385   delete [] sharingProcs;
00386 
00387   return(0);
00388 }
00389 
00390 int init_constraints(FEI* fei, HexBeam& hexcube, int& firstLocalCRID)
00391 {
00392   int numCRs = hexcube.getNumCRs();
00393   if (numCRs < 1) {
00394     return(0);
00395   }
00396 
00397   int numNodesPerCR = hexcube.getNumNodesPerCR();
00398   int* crnodes_1d = new int[numCRs*numNodesPerCR];
00399   int** crNodes = new int*[numCRs];
00400   int i, offset = 0;
00401   for(i=0; i<numCRs; ++i) {
00402     crNodes[i] = &(crnodes_1d[offset]);
00403     offset += numNodesPerCR;
00404   }
00405 
00406   CHK_ERR( hexcube.getCRNodes(crNodes) );
00407 
00408   int crID;
00409   int* fieldIDs = new int[numNodesPerCR];
00410   for(i=0; i<numNodesPerCR; ++i) fieldIDs[i] = 0;
00411 
00412   for(i=0; i<numCRs; ++i) {
00413     CHK_ERR( fei->initCRMult(numNodesPerCR, crNodes[i], fieldIDs, crID) );
00414 //     FEI_COUT << "crID: " << crID << ", nodes: ";
00415 //     for(int j=0; j<numNodesPerCR; ++j) {
00416 //       FEI_COUT << crNodes[i][j] << " ";
00417 //     }
00418 //     FEI_COUT << FEI_ENDL;
00419 
00420     if (i == 0) {
00421       firstLocalCRID = crID;
00422     }
00423   }
00424 
00425   delete [] crnodes_1d;
00426   delete [] crNodes;
00427   delete [] fieldIDs;
00428 
00429   return(0);
00430 }
00431 
00432 int load_constraints(FEI* fei, HexBeam& hexcube, int firstLocalCRID)
00433 {
00434   int numCRs = hexcube.getNumCRs();
00435   if (numCRs < 1) {
00436     return(0);
00437   }
00438 
00439   int numNodesPerCR = hexcube.getNumNodesPerCR();
00440   int* crnodes_1d = new int[numCRs*numNodesPerCR];
00441   int** crNodes = new int*[numCRs];
00442   int i, offset = 0;
00443   for(i=0; i<numCRs; ++i) {
00444     crNodes[i] = &(crnodes_1d[offset]);
00445     offset += numNodesPerCR;
00446   }
00447 
00448   CHK_ERR( hexcube.getCRNodes(crNodes) );
00449 
00450   int* fieldIDs = new int[numNodesPerCR];
00451   for(i=0; i<numNodesPerCR; ++i) fieldIDs[i] = 0;
00452 
00453   int fieldSize = hexcube.numDofPerNode();
00454   double* weights = new double[fieldSize*numNodesPerCR];
00455 
00456   for(i=0; i<fieldSize*numNodesPerCR; ++i) weights[i] = 0.0;
00457   weights[0] = -1.0;
00458   weights[fieldSize] = 1.0;
00459   double rhsValue = 0.0;
00460 
00461   for(i=0; i<numCRs; ++i) {
00462     CHK_ERR( fei->loadCRMult(firstLocalCRID+i,
00463            numNodesPerCR, crNodes[i], fieldIDs,
00464            weights, rhsValue) );
00465   }
00466 
00467   delete [] crnodes_1d;
00468   delete [] crNodes;
00469   delete [] fieldIDs;
00470   delete [] weights;
00471 
00472   return(0);
00473 }
00474 
00475 int load_elem_data(FEI* fei, HexBeam& hexcube)
00476 {
00477   int blockID = 0;
00478   int numLocalElems = hexcube.localNumElems_;
00479   int firstLocalElem = hexcube.firstLocalElem_;
00480   int nodesPerElem = hexcube.numNodesPerElem();
00481   int fieldSize = hexcube.numDofPerNode();
00482 
00483   int len = nodesPerElem*fieldSize;
00484   double* elemMat = new double[len*len];
00485   double** elemMat2D = new double*[len];
00486   double* elemVec = new double[len];
00487 
00488   int* nodeIDs = new int[nodesPerElem];
00489 
00490   if (elemMat == NULL || elemMat2D == NULL || elemVec == NULL ||
00491       nodeIDs == NULL) {
00492     return(-1);
00493   }
00494 
00495   for(int j=0; j<len; ++j) {
00496     elemMat2D[j] = &(elemMat[j*len]);
00497   }
00498 
00499   CHK_ERR( hexcube.getElemStiffnessMatrix(firstLocalElem, elemMat) );
00500   CHK_ERR( hexcube.getElemLoadVector(firstLocalElem, elemVec) );
00501 
00502   for(int i=0; i<numLocalElems; ++i) {
00503     CHK_ERR( hexcube.getElemConnectivity(firstLocalElem+i, nodeIDs) );
00504 
00505     CHK_ERR( fei->sumInElemMatrix(blockID, firstLocalElem+i,
00506           nodeIDs, elemMat2D, FEI_DENSE_ROW) );
00507 
00508     CHK_ERR( fei->sumInElemRHS(blockID, firstLocalElem+i, nodeIDs, elemVec) );
00509   }
00510 
00511   delete [] elemMat;
00512   delete [] elemMat2D;
00513   delete [] elemVec;
00514   delete [] nodeIDs;
00515 
00516   return(0);
00517 }
00518 
00519 int load_BC_data(FEI* fei, HexBeam& hexcube)
00520 {
00521   int numBCNodes = hexcube.getNumBCNodes();
00522   if (numBCNodes == 0) {
00523     return(0);
00524   }
00525 
00526   int* nodeIDs = new int[numBCNodes];
00527 
00528   int fieldID = 0;
00529 
00530   int* offsetsIntoField = new int[numBCNodes];
00531   double* prescribed_vals = new double[numBCNodes];
00532 
00533   CHK_ERR( hexcube.getBCNodes(numBCNodes, nodeIDs) );
00534 
00535   CHK_ERR( hexcube.getBCValues(numBCNodes, offsetsIntoField, prescribed_vals) );
00536 
00537   CHK_ERR( fei->loadNodeBCs(numBCNodes, nodeIDs,
00538           fieldID, offsetsIntoField, prescribed_vals) );
00539 
00540   delete [] nodeIDs;
00541   delete [] offsetsIntoField;
00542   delete [] prescribed_vals;
00543 
00544   return(0);
00545 }
00546 
00547 int init_elem_connectivities(fei::MatrixGraph* matrixGraph, HexBeam& hexcube)
00548 {
00549   int numLocalElems = hexcube.localNumElems_;
00550   int firstLocalElem = hexcube.firstLocalElem_;
00551   int nodesPerElem = hexcube.numNodesPerElem();
00552   int fieldID = 0;
00553 //  int fieldSize = hexcube.numDofPerNode();
00554   int nodeIDType = 0;
00555 
00556 
00557   int patternID = 0;
00558 //  if (fieldSize > 1) {
00559     patternID = matrixGraph->definePattern(nodesPerElem,
00560              nodeIDType, fieldID);
00561 //  }
00562 //  else {
00563 //    //if fieldSize == 1, let's not even bother associating a field with
00564 //    //our mesh-nodes. fei:: objects assume that identifiers without an
00565 //    //associated field always have exactly one degree-of-freedom.
00566 //    //
00567 //    patternID = matrixGraph->definePattern(nodesPerElem, nodeIDType);
00568 //  }
00569 
00570   int blockID = 0;
00571   CHK_ERR( matrixGraph->initConnectivityBlock(blockID, numLocalElems, patternID));
00572 
00573   int* nodeIDs = new int[nodesPerElem];
00574   if (nodeIDs == NULL) return(-1);
00575 
00576   for(int i=0; i<numLocalElems; ++i) {
00577     CHK_ERR( hexcube.getElemConnectivity(firstLocalElem+i, nodeIDs) );
00578 
00579     CHK_ERR( matrixGraph->initConnectivity(blockID, firstLocalElem+i, nodeIDs) );
00580   }
00581 
00582   delete [] nodeIDs;
00583 
00584   return(0);
00585 }
00586 
00587 int init_shared_nodes(fei::MatrixGraph* matrixGraph, HexBeam& hexcube)
00588 {
00589   int numSharedNodes = hexcube.getNumSharedNodes();
00590   if (numSharedNodes == 0) {
00591     return(0);
00592   }
00593 
00594   int* sharedNodes = NULL;
00595   int* numSharingProcsPerNode = NULL;
00596   int** sharingProcs = NULL;
00597   if (numSharedNodes > 0) {
00598     CHK_ERR( hexcube.getSharedNodes(numSharedNodes,
00599             sharedNodes, numSharingProcsPerNode,
00600             sharingProcs) );
00601   }
00602 
00603   fei::SharedPtr<fei::VectorSpace> nodeSpace = matrixGraph->getRowSpace();
00604 
00605   int nodeIDType = 0;
00606 
00607   CHK_ERR( nodeSpace->initSharedIDs(numSharedNodes, nodeIDType,
00608             sharedNodes, numSharingProcsPerNode,
00609             sharingProcs) );
00610 
00611   delete [] sharedNodes;
00612   delete [] numSharingProcsPerNode;
00613   delete [] sharingProcs[0];
00614   delete [] sharingProcs;
00615 
00616   return(0);
00617 }
00618 
00619 int init_constraints(fei::MatrixGraph* matrixGraph, HexBeam& hexcube,
00620          int localProc, int& firstLocalCRID)
00621 {
00622   int numCRs = hexcube.getNumCRs();
00623   if (numCRs < 1) {
00624     return(0);
00625   }
00626 
00627   int numNodesPerCR = hexcube.getNumNodesPerCR();
00628   int* crnodes_1d = new int[numCRs*numNodesPerCR];
00629   int** crNodes = new int*[numCRs];
00630   int i, offset = 0;
00631   for(i=0; i<numCRs; ++i) {
00632     crNodes[i] = &(crnodes_1d[offset]);
00633     offset += numNodesPerCR;
00634   }
00635 
00636   CHK_ERR( hexcube.getCRNodes(crNodes) );
00637 
00638   int crID = localProc*100000;
00639   firstLocalCRID = crID;
00640 
00641   int nodeIDType = 0;
00642   int crIDType = 1;
00643 
00644   int* fieldIDs = new int[numNodesPerCR];
00645   int* idTypes = new int[numNodesPerCR];
00646   for(i=0; i<numNodesPerCR; ++i) {
00647     fieldIDs[i] = 0;
00648     idTypes[i] = nodeIDType;
00649   }
00650 
00651   for(i=0; i<numCRs; ++i) {
00652     CHK_ERR( matrixGraph->initLagrangeConstraint(crID+i, crIDType,
00653              numNodesPerCR,
00654              idTypes, crNodes[i],
00655              fieldIDs) );
00656 //     FEI_COUT << "crID: " << crID << ", nodes: ";
00657 //     for(int j=0; j<numNodesPerCR; ++j) {
00658 //       FEI_COUT << crNodes[i][j] << " ";
00659 //     }
00660 //     FEI_COUT << FEI_ENDL;
00661 
00662   }
00663 
00664   delete [] crnodes_1d;
00665   delete [] crNodes;
00666   delete [] fieldIDs;
00667 
00668   return(0);
00669 }
00670 
00671 int init_slave_constraints(fei::MatrixGraph* matrixGraph, HexBeam& hexcube)
00672 {
00673   int numCRs = hexcube.getNumCRs();
00674   if (numCRs < 1) {
00675     return(0);
00676   }
00677 
00678   int numNodesPerCR = hexcube.getNumNodesPerCR();
00679   int* crnodes_1d = new int[numCRs*numNodesPerCR];
00680   int** crNodes = new int*[numCRs];
00681   int i, offset = 0;
00682   for(i=0; i<numCRs; ++i) {
00683     crNodes[i] = &(crnodes_1d[offset]);
00684     offset += numNodesPerCR;
00685   }
00686 
00687   CHK_ERR( hexcube.getCRNodes(crNodes) );
00688 
00689   int nodeIDType = 0;
00690 
00691   int* fieldIDs = new int[numNodesPerCR];
00692   int* idTypes = new int[numNodesPerCR];
00693   for(i=0; i<numNodesPerCR; ++i) {
00694     fieldIDs[i] = 0;
00695     idTypes[i] = nodeIDType;
00696   }
00697 
00698   int fieldSize = hexcube.numDofPerNode();
00699   double* weights = new double[fieldSize*numNodesPerCR];
00700 
00701   for(i=0; i<fieldSize*numNodesPerCR; ++i) weights[i] = 0.0;
00702   weights[0] = -1.0;
00703   weights[fieldSize] = 1.0;
00704   double rhsValue = 0.0;
00705   int offsetOfSlave = 0;
00706   int offsetIntoSlaveField = 0;
00707 
00708   for(i=0; i<numCRs; ++i) {
00709     CHK_ERR( matrixGraph->initSlaveConstraint(numNodesPerCR,
00710                 idTypes,
00711                 crNodes[i],
00712                 fieldIDs,
00713                 offsetOfSlave,
00714                 offsetIntoSlaveField,
00715                 weights,
00716                 rhsValue) );
00717   }
00718 
00719   delete [] crnodes_1d;
00720   delete [] crNodes;
00721   delete [] fieldIDs;
00722   delete [] weights;
00723 
00724   return(0);
00725 }
00726 
00727 int load_elem_data(fei::MatrixGraph* matrixGraph,
00728        fei::Matrix* mat,
00729        fei::Vector* rhs,
00730        HexBeam& hexcube)
00731 {
00732   int blockID = 0;
00733   int numLocalElems = hexcube.localNumElems_;
00734   int firstLocalElem = hexcube.firstLocalElem_;
00735   int nodesPerElem = hexcube.numNodesPerElem();
00736   int fieldSize = hexcube.numDofPerNode();
00737 
00738   int len = nodesPerElem*fieldSize;
00739   double* elemMat = new double[len*len];
00740   double** elemMat2D = new double*[len];
00741   double* elemVec = new double[len];
00742 
00743   if (elemMat == NULL || elemMat2D == NULL || elemVec == NULL) {
00744     return(-1);
00745   }
00746 
00747   for(int j=0; j<len; ++j) {
00748     elemMat2D[j] = &(elemMat[j*len]);
00749   }
00750 
00751   CHK_ERR( hexcube.getElemStiffnessMatrix(firstLocalElem, elemMat) );
00752   CHK_ERR( hexcube.getElemLoadVector(firstLocalElem, elemVec) );
00753 
00754   bool block_matrix = mat->usingBlockEntryStorage();
00755   std::vector<int> indices(len);
00756 
00757   if (block_matrix) {
00758     mat->getMatrixGraph()->setIndicesMode(fei::MatrixGraph::POINT_ENTRY_GRAPH);
00759   }
00760 
00761   for(int i=0; i<numLocalElems; ++i) {
00762     CHK_ERR( mat->getMatrixGraph()->getConnectivityIndices(blockID,
00763                                                            firstLocalElem+i,
00764                                                            len, &indices[0],
00765                                                            len) );
00766     CHK_ERR( mat->sumIn(len, &indices[0], len, &indices[0],
00767                         elemMat2D, FEI_DENSE_COL) );
00768     CHK_ERR( rhs->sumIn(len, &indices[0], elemVec, 0) );
00769   }
00770 
00771   if (block_matrix) {
00772     mat->getMatrixGraph()->setIndicesMode(fei::MatrixGraph::BLOCK_ENTRY_GRAPH);
00773   }
00774   delete [] elemMat;
00775   delete [] elemMat2D;
00776   delete [] elemVec;
00777 
00778   return(0);
00779 }
00780 
00781 int load_constraints(fei::LinearSystem* linSys, HexBeam& hexcube,
00782          int firstLocalCRID)
00783 {
00784   int numCRs = hexcube.getNumCRs();
00785   if (numCRs < 1) {
00786     return(0);
00787   }
00788 
00789   int numNodesPerCR = hexcube.getNumNodesPerCR();
00790 
00791   int fieldSize = hexcube.numDofPerNode();
00792   double* weights = new double[fieldSize*numNodesPerCR];
00793 
00794   int i;
00795   for(i=0; i<fieldSize*numNodesPerCR; ++i) weights[i] = 0.0;
00796   weights[0] = -1.0;
00797   weights[fieldSize] = 1.0;
00798   double rhsValue = 0.0;
00799 
00800   for(i=0; i<numCRs; ++i) {
00801     CHK_ERR( linSys->loadLagrangeConstraint(firstLocalCRID+i,
00802               weights, rhsValue) );
00803   }
00804 
00805   delete [] weights;
00806 
00807   return(0);
00808 }
00809 
00810 int load_BC_data(fei::LinearSystem* linSys, HexBeam& hexcube)
00811 {
00812   int numBCNodes = hexcube.getNumBCNodes();
00813   if (numBCNodes == 0) {
00814     return(0);
00815   }
00816 
00817   int* nodeIDs = new int[numBCNodes];
00818 
00819   int fieldID = 0;
00820   int nodeIDType = 0;
00821 
00822   int* offsetsIntoField = new int[numBCNodes];
00823   double* prescribed_vals = new double[numBCNodes];
00824 
00825   CHK_ERR( hexcube.getBCNodes(numBCNodes, nodeIDs) );
00826 
00827   CHK_ERR( hexcube.getBCValues(numBCNodes, offsetsIntoField, prescribed_vals) );
00828 
00829   CHK_ERR( linSys->loadEssentialBCs(numBCNodes, nodeIDs,
00830             nodeIDType, fieldID,
00831             offsetsIntoField, prescribed_vals) );
00832 
00833   delete [] offsetsIntoField;
00834   delete [] prescribed_vals;
00835   delete [] nodeIDs;
00836 
00837   return(0);
00838 }
00839 
00840 }//namespace HexBeam_Functions
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends