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