PoissonData.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 //
00010 //This is a class that will exercise FEI implementations.
00011 //
00012 
00013 #include <fei_macros.hpp>
00014 #include <fei_defs.h>
00015 #include <fei_CSVec.hpp>
00016 
00017 #include <test_utils/Poisson_Elem.hpp>
00018 
00019 #include <test_utils/PoissonData.hpp>
00020 
00021 #include <cstdlib>
00022 #include <cmath>
00023 
00024 static int int_sqrt(int x) {
00025 //a not-safe function for taking the sqrt of a square int.
00026     return((int)std::ceil(std::sqrt((double)x)));
00027 }
00028 
00029 //==============================================================================
00030 PoissonData::PoissonData(int L,
00031                        int numProcs, int localProc, int outputLevel)
00032 {
00033   //
00034   //PoissonData constructor.
00035   //
00036   //Arguments:
00037   //
00038   // L:                global square size (number-of-elements along side)
00039   // numProcs:         number of processors participating in this FEI test.
00040   // localProc:        local processor number.
00041   // outputLevel:      affects the amount of screen output.
00042   //
00043 
00044   L_ = L;
00045     
00046   startElement_ = 0;
00047   numLocalElements_ = 0;
00048 
00049   numProcs_ = numProcs;
00050   localProc_ = localProc;
00051   outputLevel_ = outputLevel;
00052 
00053   check1();
00054 
00055   elem_ = new Poisson_Elem();
00056   int err = elem_->allocateInternals(1);
00057   err += elem_->allocateLoad(1);
00058   err += elem_->allocateStiffness(1);
00059   if (err) messageAbort("Allocation error in element.");
00060 
00061   fieldArraysAllocated_ = false;
00062   elemIDsAllocated_ = false;
00063 
00064   numFields_ = NULL;
00065   fieldIDs_ = NULL;
00066 
00067   elemIDs_ = NULL;
00068 
00069   calculateDistribution();
00070 
00071   numElemBlocks_ = 1;
00072   elemBlockID_ = (GlobalID)0;
00073   elemSetID_ = 0;
00074   elemFormat_ = 0;
00075 
00076   nodesPerElement_ = 4;
00077   fieldsPerNode_ = 1;
00078 
00079   initializeFieldStuff();
00080 }
00081 
00082 //==============================================================================
00083 PoissonData::~PoissonData() {
00084 //
00085 //Destructor -- delete any allocated memory.
00086 //
00087 
00088     deleteFieldArrays();
00089 
00090     if (elemIDsAllocated_) delete [] elemIDs_;
00091 
00092     elem_->deleteMemory();
00093     delete elem_;
00094 }
00095 
00096 //==============================================================================
00097 void PoissonData::check1() {
00098 //
00099 //Private function to be called from the constructor, simply makes sure that
00100 //the constructor's input arguments were reasonable.
00101 //
00102 //If they aren't, a message is printed on standard err, and abort() is called.
00103 //
00104     if (L_ <= 0)                 messageAbort("bar length L <= 0.");
00105     if (numProcs_ <= 0)          messageAbort("numProcs <= 0.");
00106     if (L_%int_sqrt(numProcs_)) 
00107         messageAbort("L must be an integer multiple of sqrt(numProcs).");
00108     if (localProc_ < 0)          messageAbort("localProc < 0.");
00109     if (localProc_ >= numProcs_) messageAbort("localProc >= numProcs.");
00110     if (outputLevel_ < 0)        messageAbort("outputLevel < 0.");
00111 }
00112 
00113 //==============================================================================
00114 void PoissonData::calculateDistribution() {
00115 //
00116 //Calculate which elements this processor owns. The element domain is a
00117 //square, and we can assume that sqrt(numProcs_) divides evenly into
00118 //L_. We're working with a (logically) 2D processor arrangement.
00119 //Furthermore, the logical processor layout is such that processor 0 is at
00120 //the bottom left corner of a 2D grid, and a side of the grid is of length
00121 //sqrt(numProcs_). The element domain is numbered such that element 1 is at
00122 //the bottom left corner of the square, and element numbers increase from
00123 //left to right. i.e., element 1 is in position (1,1), element L is in
00124 //position (1,L), element L+1 is in position (2,1).
00125 //
00126 //Use 1-based numbering for the elements and the x- and y- coordinates in
00127 //the element grid, but use 0-based numbering for processor IDs and the
00128 //coordinates in the processor grid.
00129 //
00130     numLocalElements_ = (L_*L_)/numProcs_;
00131 
00132     elemIDs_ = new GlobalID[numLocalElements_];
00133     if (!elemIDs_) messageAbort("ERROR allocating elemIDs_.");
00134     elemIDsAllocated_ = true;
00135 
00136     //0-based x-coordinate of this processor in the 2D processor grid.
00137     procX_ = localProc_%int_sqrt(numProcs_);
00138 
00139     //0-based maximum processor x-coordinate.
00140     maxProcX_ = int_sqrt(numProcs_) - 1;
00141 
00142     //0-based y-coordinate of this processor in the 2D processor grid.
00143     procY_ = localProc_/int_sqrt(numProcs_);
00144 
00145     //0-based maximum processor y-coordinate.
00146     maxProcY_ = int_sqrt(numProcs_) - 1;
00147 
00148     int sqrtElems = int_sqrt(numLocalElements_);
00149     int sqrtProcs = int_sqrt(numProcs_);
00150 
00151     //1-based first-element-on-this-processor
00152     startElement_ = 1 + procY_*sqrtProcs*numLocalElements_ + procX_*sqrtElems;
00153 
00154     if (outputLevel_>1) {
00155         FEI_COUT << localProc_ << ", calcDist.: numLocalElements: " 
00156              << numLocalElements_ << ", startElement: " << startElement_ 
00157              << FEI_ENDL;
00158         FEI_COUT << localProc_ << ", procX: " << procX_ << ", procY_: " << procY_
00159              << ", maxProcX: " << maxProcX_ << ", maxProcY: " << maxProcY_
00160              << FEI_ENDL;
00161     }
00162 
00163     int offset = 0;
00164     for(int i=0; i<sqrtElems; i++) {
00165         for(int j=0; j<sqrtElems; j++) {
00166             elemIDs_[offset] = (GlobalID)(startElement_ + i*L_ + j);
00167             offset++;
00168         }
00169     }
00170 }
00171 
00172 //==============================================================================
00173 void PoissonData::messageAbort(const char* message) {
00174     FEI_CERR << FEI_ENDL << "PoissonData: " << message 
00175          << FEI_ENDL << "  Aborting." << FEI_ENDL;
00176     std::abort();
00177 }
00178 
00179 //==============================================================================
00180 GlobalID* PoissonData::getElementConnectivity(GlobalID elemID)
00181 {
00182    //set the elemID on the internal Poisson_Elem instance.
00183    elem_->setElemID(elemID);
00184    elem_->setElemLength(1.0/L_);
00185    elem_->setTotalLength(1.0);
00186 
00187    //now get a pointer to the element's connectivity array and
00188    //calculate that connectivity (in place).
00189    int size = 0;
00190    GlobalID* elemConn = elem_->getElemConnPtr(size);
00191    if (size == 0) messageAbort("loadElements: bad conn ptr.");
00192 
00193    calculateConnectivity(elemConn, size, elemID);
00194 
00195    return(elemConn);
00196 }
00197 
00198 //==============================================================================
00199 double** PoissonData::getElemStiffness(GlobalID elemID)
00200 {
00201    elem_->setElemID(elemID);
00202    elem_->setElemLength(1.0/L_);
00203    elem_->setTotalLength(1.0);
00204 
00205    //now get a pointer to this element's connectivity array and
00206    //calculate that connectivity (in place).
00207    int size = 0;
00208    GlobalID* elemConn = elem_->getElemConnPtr(size);
00209    if (size == 0) messageAbort("loadElemStiffnesses: bad conn ptr.");
00210 
00211    calculateConnectivity(elemConn, size, elemID);
00212 
00213    elem_->calculateCoords();
00214 
00215    if (outputLevel_>1) {
00216       double* x = elem_->getNodalX(size);
00217       double* y = elem_->getNodalY(size);
00218       FEI_COUT << localProc_ << ", elemID " << elemID << ", nodes: ";
00219       for(int j=0; j<size; j++) {
00220          FEI_COUT << elemConn[j] << " ";
00221          FEI_COUT << "("<<x[j]<<","<<y[j]<<") ";
00222       }
00223       FEI_COUT << FEI_ENDL;
00224    }
00225 
00226    elem_->calculateStiffness();
00227 
00228    return( elem_->getElemStiff(size) );
00229 }
00230 
00231 //==============================================================================
00232 double* PoissonData::getElemLoad(GlobalID elemID)
00233 {
00234    elem_->setElemID(elemID);
00235    elem_->setElemLength(1.0/L_);
00236    elem_->setTotalLength(1.0);
00237 
00238    //now get a pointer to this element's connectivity array and
00239    //calculate that connectivity (in place).
00240    int size = 0;
00241    GlobalID* elemConn = elem_->getElemConnPtr(size);
00242    if (size == 0) messageAbort("loadElemLoads: bad conn ptr.");
00243 
00244    calculateConnectivity(elemConn, size, elemID);
00245 
00246    elem_->calculateCoords();
00247 
00248    if (outputLevel_>1) {
00249       double* x = elem_->getNodalX(size);
00250       double* y = elem_->getNodalY(size);
00251       FEI_COUT << localProc_ << ", elemID " << elemID << ", nodes: ";
00252       for(int j=0; j<size; j++) {
00253          FEI_COUT << elemConn[j] << " ";
00254          FEI_COUT << "("<<x[j]<<","<<y[j]<<") ";
00255       }
00256       FEI_COUT << FEI_ENDL;
00257    }
00258 
00259    elem_->calculateLoad();
00260 
00261    return( elem_->getElemLoad(size));
00262 
00263 }
00264 
00265 //==============================================================================
00266 void PoissonData::calculateConnectivity(GlobalID* conn, int size,
00267                                         GlobalID elemID) {
00268 //
00269 //Calculate a single element's connectivity array -- the list of nodes
00270 //that it 'contains'.
00271 //
00272 //Note that we're assuming the element is a 2D square.
00273 //
00274     //elemX will be the global 'x-coordinate' of this element in the square. The
00275     //'origin' is the lower-left corner of the bar, which is element 1,
00276     //and it is in position 1,1.
00277     int elemX = (int)elemID%L_;
00278     if (elemX == 0) elemX = L_;
00279 
00280     //elemY will be the global (1-based) 'y-coordinate'.
00281     int elemY = ((int)elemID - elemX)/L_ + 1;
00282 
00283     //These are the four nodes for this element.
00284     GlobalID lowerLeft = elemID + (GlobalID)(elemY-1);
00285     GlobalID lowerRight = lowerLeft + (GlobalID)1;
00286     GlobalID upperRight = lowerRight + (GlobalID)(L_+1);
00287     GlobalID upperLeft = upperRight - (GlobalID)1;
00288 
00289     (void)size;
00290 
00291     //now fill the connectivity array. We'll always fill the connectivity
00292     //array with the lower left node first, and then proceed counter-clockwise.
00293     conn[0] = lowerLeft;
00294     conn[1] = lowerRight;
00295     conn[2] = upperRight;
00296     conn[3] = upperLeft;
00297 }
00298 
00299 //==============================================================================
00300 void PoissonData::initializeFieldStuff() {
00301 //
00302 //Set up the field-descriptor variables that will be passed
00303 //to the FEI's initFields function, beginInitElemBlock function, etc.
00304 //
00305 //Note we've hardwired 1 dof per field.
00306 //
00307     fieldSize_ = 1;
00308     numFields_ = new int[nodesPerElement_];
00309     fieldIDs_ = new int*[nodesPerElement_];
00310     for(int i=0; i<nodesPerElement_; i++) {
00311         numFields_[i] = fieldsPerNode_;
00312         fieldIDs_[i] = new int[fieldsPerNode_];
00313         for(int j=0; j<fieldsPerNode_; j++) {
00314             fieldIDs_[i][j] = fieldsPerNode_;
00315         }
00316     }
00317     fieldArraysAllocated_ = true;
00318 }
00319 
00320 //==============================================================================
00321 void PoissonData::deleteFieldArrays() {
00322 
00323     if (fieldArraysAllocated_) {
00324 
00325         for(int i=0; i<nodesPerElement_; i++) {
00326             delete [] fieldIDs_[i];
00327         }
00328 
00329         delete [] fieldIDs_;
00330         delete [] numFields_;
00331     }
00332     fieldArraysAllocated_ = false;
00333 }
00334 
00335 //==============================================================================
00336 void PoissonData::getLeftSharedNodes(int& numShared, GlobalID* sharedNodeIDs,
00337                                      int* numProcsPerSharedNode,
00338                                      int** sharingProcs) {
00339 //
00340 //This function decides whether any of the nodes along the left edge,
00341 //including the top node but not the bottom node, are shared. It also
00342 //decides which processors the nodes are shared with.
00343 //
00344 
00345     if (numProcs_ == 1) {
00346         numShared = 0;
00347         return;
00348     }
00349 
00350     if (procX_ == 0) {
00351         //if this proc is on the left edge of the square...
00352 
00353         if (procY_ < maxProcY_) {
00354             //if this proc is not the top left proc...
00355 
00356             numShared = 1;
00357 
00358             int topLeftElemIndex = numLocalElements_ -
00359                                int_sqrt(numLocalElements_);
00360 
00361             elem_->setElemID(elemIDs_[topLeftElemIndex]);
00362 
00363             //now get a pointer to this element's connectivity array and
00364             //calculate that connectivity (in place).
00365             int size = 0;
00366             GlobalID* elemConn = elem_->getElemConnPtr(size);
00367             if (size == 0) messageAbort("loadElements: bad conn ptr.");
00368 
00369             calculateConnectivity(elemConn, size, elemIDs_[topLeftElemIndex]);
00370 
00371             sharedNodeIDs[0] = elemConn[3]; //elem's top left node is node 3
00372             numProcsPerSharedNode[0] = 2;
00373             sharingProcs[0][0] = localProc_;
00374             sharingProcs[0][1] = localProc_ + int_sqrt(numProcs_);
00375 
00376             return;
00377         }
00378         else {
00379             //else this proc is the top left proc...
00380             numShared = 0;
00381         }
00382     }
00383     else {
00384         //else this proc is not on the left edge of the square...
00385 
00386         numShared = int_sqrt(numLocalElements_);
00387         int lowerLeftElemIndex = 0;
00388 
00389         int sqrtElems = int_sqrt(numLocalElements_);
00390 
00391         int shOffset = 0;
00392         for(int i=0; i<sqrtElems; i++){
00393             //stride up the left edge of the local elements...
00394             int size=0;
00395 
00396             int elemIndex = lowerLeftElemIndex+i*sqrtElems;
00397 
00398             elem_->setElemID(elemIDs_[elemIndex]);
00399 
00400             //now get a pointer to this element's connectivity array and
00401             //calculate that connectivity (in place).
00402             GlobalID* nodes = elem_->getElemConnPtr(size);
00403             if (size == 0) messageAbort(": bad conn ptr.");
00404       
00405             calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
00406 
00407             //now put in the top left node
00408             sharedNodeIDs[shOffset] = nodes[3];
00409             sharingProcs[shOffset][0] = localProc_-1;
00410             sharingProcs[shOffset][1] = localProc_;
00411             numProcsPerSharedNode[shOffset++] = 2;
00412         }
00413 
00414         if (procY_ < maxProcY_) {
00415             //if this proc isn't on the top edge, the upper left node (the
00416             //last one we put into the shared node list) is shared by 4 procs.
00417             shOffset--;
00418             numProcsPerSharedNode[shOffset] = 4;
00419             sharingProcs[shOffset][2] = localProc_ + int_sqrt(numProcs_);
00420             sharingProcs[shOffset][3] = sharingProcs[shOffset][2] - 1;
00421         }
00422     }
00423 }
00424 
00425 //==============================================================================
00426 void PoissonData::getRightSharedNodes(int& numShared, GlobalID* sharedNodeIDs,
00427                                      int* numProcsPerSharedNode,
00428                                      int** sharingProcs) {
00429 //
00430 //This function decides whether any of the nodes along the right edge,
00431 //including the bottom node but not the top node, are shared. It also
00432 //decides which processors the nodes are shared with.
00433 //
00434 
00435     if (numProcs_ == 1) {
00436         numShared = 0;
00437         return;
00438     }
00439 
00440     if (procX_ == maxProcX_) {
00441         //if this proc is on the right edge of the square...
00442 
00443         if (procY_ > 0) {
00444             //if this proc is not the bottom right proc...
00445 
00446             numShared = 1;
00447 
00448             int lowerRightElemIndex = int_sqrt(numLocalElements_) - 1;
00449 
00450             elem_->setElemID(elemIDs_[lowerRightElemIndex]);
00451 
00452             //now get a pointer to this element's connectivity array and
00453             //calculate that connectivity (in place).
00454             int size;
00455             GlobalID* nodes = elem_->getElemConnPtr(size);
00456             if (size == 0) messageAbort(": bad conn ptr.");
00457       
00458             calculateConnectivity(nodes, size, elemIDs_[lowerRightElemIndex]);
00459 
00460             sharedNodeIDs[0] = nodes[1]; //elem's bottom right node is node 1
00461             numProcsPerSharedNode[0] = 2;
00462             sharingProcs[0][0] = localProc_;
00463             sharingProcs[0][1] = localProc_ - int_sqrt(numProcs_);
00464 
00465             return;
00466         }
00467         else {
00468             //else this proc is the bottom right proc...
00469             numShared = 0;
00470         }
00471     }
00472     else {
00473         //else this proc is not on the right edge of the square...
00474 
00475         numShared = int_sqrt(numLocalElements_);
00476         int upperRightElemIndex = numLocalElements_ - 1;
00477 
00478         int sqrtElems = int_sqrt(numLocalElements_);
00479 
00480         int shOffset = 0;
00481         for(int i=0; i<sqrtElems; i++){
00482             //stride down the right edge of the local elements...
00483             int size=0;
00484             int elemIndex = upperRightElemIndex-i*sqrtElems;
00485             elem_->setElemID(elemIDs_[elemIndex]);
00486 
00487             //now get a pointer to this element's connectivity array and
00488             //calculate that connectivity (in place).
00489             GlobalID* nodes = elem_->getElemConnPtr(size);
00490             if (size == 0) messageAbort(": bad conn ptr.");
00491 
00492             calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
00493 
00494             //now put in the lower right node
00495             sharedNodeIDs[shOffset] = nodes[1];
00496             sharingProcs[shOffset][0] = localProc_+1;
00497             sharingProcs[shOffset][1] = localProc_;
00498             numProcsPerSharedNode[shOffset++] = 2;
00499         }
00500 
00501         if (procY_ > 0) {
00502             //if this proc isn't on the bottom edge, the lower right node (the
00503             //last one we put into the shared node list) is shared by 4 procs.
00504             shOffset--;
00505             numProcsPerSharedNode[shOffset] = 4;
00506             sharingProcs[shOffset][2] = localProc_ - int_sqrt(numProcs_);
00507             sharingProcs[shOffset][3] = sharingProcs[shOffset][2] + 1;
00508         }
00509     }
00510 }
00511 
00512 //==============================================================================
00513 void PoissonData::getTopSharedNodes(int& numShared, GlobalID* sharedNodeIDs,
00514                                      int* numProcsPerSharedNode,
00515                                      int** sharingProcs) {
00516 //
00517 //This function decides whether any of the nodes along the top edge,
00518 //including the right node but not the left node, are shared. It also
00519 //decides which processors the nodes are shared with.
00520 //
00521 
00522     if (numProcs_ == 1) {
00523         numShared = 0;
00524         return;
00525     }
00526 
00527     if (procY_ == maxProcY_) {
00528         //if this proc is on the top edge of the square...
00529 
00530         if (procX_ < maxProcX_) {
00531             //if this proc is not the top right proc...
00532 
00533             numShared = 1;
00534 
00535             int elemIndex = numLocalElements_ - 1;
00536 
00537             elem_->setElemID(elemIDs_[elemIndex]);
00538 
00539             //now get a pointer to this element's connectivity array and
00540             //calculate that connectivity (in place).
00541             int size;
00542             GlobalID* nodes = elem_->getElemConnPtr(size);
00543             if (size == 0) messageAbort(": bad conn ptr.");
00544 
00545             calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
00546 
00547             sharedNodeIDs[0] = nodes[2]; //elem's top right node is node 2
00548             numProcsPerSharedNode[0] = 2;
00549             sharingProcs[0][0] = localProc_;
00550             sharingProcs[0][1] = localProc_ + 1;
00551 
00552             return;
00553         }
00554         else {
00555             //else this proc is the top right proc...
00556             numShared = 0;
00557         }
00558     }
00559     else {
00560         //else this proc is not on the top edge of the square...
00561 
00562         numShared = int_sqrt(numLocalElements_);
00563         int topLeftElemIndex = numLocalElements_ - int_sqrt(numLocalElements_);
00564 
00565         int sqrtElems = int_sqrt(numLocalElements_);
00566 
00567         int shOffset = 0;
00568         for(int i=0; i<sqrtElems; i++){
00569             //stride across the top edge of the local elements...
00570             int size=0;
00571             int elemIndex = topLeftElemIndex+i;
00572 
00573             elem_->setElemID(elemIDs_[elemIndex]);
00574 
00575             //now get a pointer to this element's connectivity array and
00576             //calculate that connectivity (in place).
00577             GlobalID* nodes = elem_->getElemConnPtr(size);
00578             if (size == 0) messageAbort(": bad conn ptr.");
00579 
00580             calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
00581 
00582             //now put in the upper right node
00583             sharedNodeIDs[shOffset] = nodes[2];
00584             sharingProcs[shOffset][0] = localProc_+int_sqrt(numProcs_);
00585             sharingProcs[shOffset][1] = localProc_;
00586             numProcsPerSharedNode[shOffset++] = 2;
00587         }
00588         if (procX_ < maxProcX_) {
00589             //if this proc isn't on the right edge, the top right node (the
00590             //last one we put into the shared node list) is shared by 4 procs.
00591             shOffset--;
00592             numProcsPerSharedNode[shOffset] = 4;
00593             sharingProcs[shOffset][2] = localProc_ + 1;
00594             sharingProcs[shOffset][3] = sharingProcs[shOffset][0] + 1;
00595         }
00596     }
00597 }
00598 
00599 //==============================================================================
00600 void PoissonData::getBottomSharedNodes(int& numShared, GlobalID* sharedNodeIDs,
00601                                      int* numProcsPerSharedNode,
00602                                      int** sharingProcs) {
00603 //
00604 //This function decides whether any of the nodes along the bottom edge,
00605 //including the left node but not the right node, are shared. It also
00606 //decides which processors the nodes are shared with.
00607 //
00608 
00609     if (numProcs_ == 1) {
00610         numShared = 0;
00611         return;
00612     }
00613 
00614     if (procY_ == 0) {
00615         //if this proc is on the bottom edge of the square...
00616 
00617         if (procX_ > 0) {
00618             //if this proc is not the bottom left proc...
00619 
00620             numShared = 1;
00621 
00622             int elemIndex = 0;
00623 
00624             elem_->setElemID(elemIDs_[elemIndex]);
00625 
00626             //now get a pointer to this element's connectivity array and
00627             //calculate that connectivity (in place).
00628             int size;
00629             GlobalID* nodes = elem_->getElemConnPtr(size);
00630             if (size == 0) messageAbort(": bad conn ptr.");
00631 
00632             calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
00633 
00634             sharedNodeIDs[0] = nodes[0]; //elem's bottom left node is node 0
00635             numProcsPerSharedNode[0] = 2;
00636             sharingProcs[0][0] = localProc_;
00637             sharingProcs[0][1] = localProc_ - 1;
00638 
00639             return;
00640         }
00641         else {
00642             //else this proc is the top right proc...
00643             numShared = 0;
00644         }
00645     }
00646     else {
00647         //else this proc is not on the bottom edge of the square...
00648 
00649         numShared = int_sqrt(numLocalElements_);
00650         int lowerRightElemIndex = int_sqrt(numLocalElements_) - 1;
00651 
00652         int sqrtElems = int_sqrt(numLocalElements_);
00653 
00654         int shOffset = 0;
00655         for(int i=0; i<sqrtElems; i++){
00656             //stride across the bottom edge of the local elements, from 
00657             //right to left...
00658             int size=0;
00659             int elemIndex = lowerRightElemIndex-i;
00660 
00661             elem_->setElemID(elemIDs_[elemIndex]);
00662 
00663             //now get a pointer to this element's connectivity array and
00664             //calculate that connectivity (in place).
00665             GlobalID* nodes = elem_->getElemConnPtr(size);
00666             if (size == 0) messageAbort(": bad conn ptr.");
00667 
00668             calculateConnectivity(nodes, size, elemIDs_[elemIndex]);
00669 
00670             //now put in the lower left node
00671             sharedNodeIDs[shOffset] = nodes[0];
00672             sharingProcs[shOffset][0] = localProc_ - int_sqrt(numProcs_);
00673             sharingProcs[shOffset][1] = localProc_;
00674             numProcsPerSharedNode[shOffset++] = 2;
00675         }
00676         if (procX_ > 0) {
00677             //if this proc isn't on the left edge, the lower left node (the
00678             //last one we put into the shared node list) is shared by 4 procs.
00679             shOffset--;
00680             numProcsPerSharedNode[shOffset] = 4;
00681             sharingProcs[shOffset][2] = localProc_ - 1;
00682             sharingProcs[shOffset][3] = sharingProcs[shOffset][0] - 1;
00683         }
00684     }
00685 }
00686 
00687 //==============================================================================
00688 void PoissonData::printSharedNodes(const char* str,
00689            int numShared, GlobalID* nodeIDs,
00690                                    int** shareProcs, int* numShareProcs)
00691 {
00692   for(int i=0; i<numShared; i++) {
00693     FEI_COUT << localProc_ << ", " << str << " node: " << (int) nodeIDs[i];
00694     FEI_COUT << ", procs: ";
00695     for(int j=0; j<numShareProcs[i]; j++) {
00696       FEI_COUT << shareProcs[i][j] << " ";
00697     }
00698     FEI_COUT << FEI_ENDL;
00699   }
00700 }
00701 
00702 //==============================================================================
00703 void PoissonData::calculateBCs() {
00704 //
00705 //This function figures out which nodes lie on the boundary. The ones that
00706 //do are added to the BC set, along with appropriate alpha/beta/gamma values.
00707 //
00708     for(int i=0; i<numLocalElements_; i++) {
00709        elem_->setElemID(elemIDs_[i]);
00710        elem_->setElemLength(1.0/L_);
00711        elem_->setTotalLength(1.0);
00712 
00713        //now get a pointer to this element's connectivity array and
00714        //calculate that connectivity (in place).
00715        int size = 0;
00716        GlobalID* nodeIDs = elem_->getElemConnPtr(size);
00717        if (size == 0) messageAbort("loadElements: bad conn ptr.");
00718 
00719        calculateConnectivity(nodeIDs, size, elemIDs_[i]);
00720 
00721        elem_->calculateCoords();
00722 
00723        double* xcoord = elem_->getNodalX(size);
00724        double* ycoord = elem_->getNodalY(size);
00725 
00726        //now loop over the nodes and see if any are on a boundary.
00727        for(int j=0; j<size; j++) {
00728           if ((std::abs(xcoord[j]) < 1.e-49) || (std::abs(xcoord[j] - 1.0) < 1.e-49) ||
00729              (std::abs(ycoord[j]) < 1.e-49) || (std::abs(ycoord[j] - 1.0) < 1.e-49)) {
00730 
00731              addBCNode(nodeIDs[j], xcoord[j], ycoord[j]);
00732           }
00733        }
00734     }
00735 }
00736 
00737 //==============================================================================
00738 void PoissonData::addBCNode(GlobalID nodeID, double x, double y){
00739 
00740   std::vector<GlobalID>::iterator
00741     iter = std::lower_bound(BCNodeIDs_.begin(), BCNodeIDs_.end(), nodeID);
00742 
00743   if (iter == BCNodeIDs_.end() || *iter != nodeID) {
00744     unsigned offset = iter - BCNodeIDs_.begin();
00745     BCNodeIDs_.insert(iter, nodeID);
00746 
00747     double bcValue = std::pow(x, 2.0) + std::pow(y, 2.0);
00748 
00749     BCValues_.insert(BCValues_.begin()+offset, bcValue);
00750   }
00751 }
00752 
00753 //==============================================================================
00754 int init_elem_connectivities(FEI* fei, PoissonData& poissonData)
00755 {
00756   //first load the information that defines this element block, and
00757   //the topology of each element in this element block.
00758 
00759   GlobalID elemBlockID = poissonData.getElemBlockID();
00760   int numLocalElements = poissonData.getNumLocalElements();
00761   int numNodesPerElement = poissonData.getNumNodesPerElement();
00762   int* numFieldsPerNode = poissonData.getNumFieldsPerNodeList();
00763   int** fieldIDsTable = poissonData.getNodalFieldIDsTable();
00764 
00765   CHK_ERR( fei->initElemBlock(elemBlockID,
00766             numLocalElements,
00767             numNodesPerElement,
00768             numFieldsPerNode,
00769             fieldIDsTable,
00770             0, // no element-centered degrees-of-freedom
00771             NULL, //null list of elem-dof fieldIDs
00772             FEI_NODE_MAJOR) );
00773 
00774   //now let's loop over all of the local elements, giving their 
00775   //nodal connectivity lists to the FEI.
00776 
00777   GlobalID* elemIDs = poissonData.getLocalElementIDs();
00778 
00779   for(int elem=0; elem<numLocalElements; elem++) {
00780     GlobalID* elemConnectivity =
00781       poissonData.getElementConnectivity(elemIDs[elem]);
00782 
00783     CHK_ERR( fei->initElem(elemBlockID, elemIDs[elem], elemConnectivity) );
00784   }
00785 
00786   return(0);
00787 }
00788 
00789 //==============================================================================
00790 int set_shared_nodes(FEI* fei, PoissonData& poissonData)
00791 {
00792    int numLocalElements = poissonData.getNumLocalElements();
00793    int maxNumSharedNodes = (int)std::sqrt((double)numLocalElements);
00794    GlobalID* sharedNodeIDs = new GlobalID[maxNumSharedNodes];
00795    int* numProcsPerSharedNode = new int[maxNumSharedNodes];
00796    int** sharingProcs = new int*[maxNumSharedNodes];
00797    for(int i=0; i<maxNumSharedNodes; i++) sharingProcs[i] = new int[4];
00798 
00799    int numShared;
00800 
00801    //first, get the shared-node data for the left edge of the local block
00802 
00803    poissonData.getLeftSharedNodes(numShared, sharedNodeIDs,
00804                                   numProcsPerSharedNode, sharingProcs);
00805 
00806    CHK_ERR( fei->initSharedNodes(numShared, sharedNodeIDs,
00807                                 numProcsPerSharedNode, sharingProcs));
00808 
00809    //now, get the shared-node data for the right edge of the local block
00810 
00811    poissonData.getRightSharedNodes(numShared, sharedNodeIDs,
00812                                   numProcsPerSharedNode, sharingProcs);
00813 
00814    CHK_ERR( fei->initSharedNodes(numShared, sharedNodeIDs,
00815                                 numProcsPerSharedNode, sharingProcs));
00816 
00817    //now, get the shared-node data for the bottom edge of the local block
00818 
00819    poissonData.getBottomSharedNodes(numShared, sharedNodeIDs,
00820                                   numProcsPerSharedNode, sharingProcs);
00821 
00822    CHK_ERR( fei->initSharedNodes(numShared, sharedNodeIDs,
00823                                 numProcsPerSharedNode, sharingProcs));
00824 
00825    //finally, get the shared-node data for the top edge of the local block
00826 
00827    poissonData.getTopSharedNodes(numShared, sharedNodeIDs,
00828                                   numProcsPerSharedNode, sharingProcs);
00829 
00830    CHK_ERR( fei->initSharedNodes(numShared, sharedNodeIDs,
00831                                 numProcsPerSharedNode, sharingProcs));
00832 
00833    for(int j=0; j<maxNumSharedNodes; j++) delete [] sharingProcs[j];
00834    delete [] sharingProcs;
00835    delete [] numProcsPerSharedNode;
00836    delete [] sharedNodeIDs;
00837 
00838    return(0);
00839 }
00840 
00841 //==============================================================================
00842 int load_elem_data(FEI* fei, PoissonData& poissonData)
00843 {
00844   GlobalID elemBlockID = poissonData.getElemBlockID();
00845   int numLocalElements = poissonData.getNumLocalElements();
00846   GlobalID* elemIDs = poissonData.getLocalElementIDs();
00847 
00848   for(int elem=0; elem<numLocalElements; elem++) {
00849     GlobalID* elemConnectivity =
00850       poissonData.getElementConnectivity(elemIDs[elem]);
00851     double** elemStiffness = poissonData.getElemStiffness(elemIDs[elem]);
00852 
00853     CHK_ERR( fei->sumInElemMatrix(elemBlockID, elemIDs[elem],
00854           elemConnectivity, elemStiffness,
00855           poissonData.getElemFormat()));
00856 
00857     double* elemLoad = poissonData.getElemLoad(elemIDs[elem]);
00858 
00859     CHK_ERR( fei->sumInElemRHS(elemBlockID, elemIDs[elem],
00860              elemConnectivity, elemLoad));
00861   }
00862 
00863   return(0);
00864 }
00865 
00866 //==============================================================================
00867 int load_elem_data_putrhs(FEI* fei, PoissonData& poissonData)
00868 {
00869   GlobalID elemBlockID = poissonData.getElemBlockID();
00870   int numLocalElements = poissonData.getNumLocalElements();
00871   GlobalID* elemIDs = poissonData.getLocalElementIDs();
00872 
00873   int numIDs = poissonData.getNumNodesPerElement();
00874 
00875   int* fieldID = poissonData.getFieldIDs();
00876 
00877   fei::CSVec rhs;
00878 
00879   for(int elem=0; elem<numLocalElements; elem++) {
00880     GlobalID* elemConnectivity =
00881       poissonData.getElementConnectivity(elemIDs[elem]);
00882     double** elemStiffness = poissonData.getElemStiffness(elemIDs[elem]);
00883 
00884     CHK_ERR( fei->sumInElemMatrix(elemBlockID, elemIDs[elem],
00885                                   elemConnectivity, elemStiffness,
00886                                   poissonData.getElemFormat()));
00887 
00888     double* elemLoad = poissonData.getElemLoad(elemIDs[elem]);
00889 
00890     for(int i=0; i<numIDs; ++i) {
00891       fei::add_entry(rhs, elemConnectivity[i], elemLoad[i]);
00892     }
00893   }
00894 
00895   fei->loadComplete();
00896 
00897   fei->putIntoRHS(0, *fieldID, rhs.size(),
00898                   &(rhs.indices()[0]), &(rhs.coefs()[0]));
00899 
00900   return(0);
00901 }
00902 
00903 //==============================================================================
00904 int load_BC_data(FEI* fei, PoissonData& poissonData)
00905 {
00906   //first, have the data object generate the BC data
00907   poissonData.calculateBCs();
00908 
00909   int numBCNodes = poissonData.getNumBCNodes();
00910   GlobalID* nodeIDs = poissonData.getBCNodeIDs();
00911   int fieldID = poissonData.getBCFieldID();
00912   double* values = poissonData.getBCValues();
00913 
00914   std::vector<int> offsets(numBCNodes, 0);
00915 
00916   CHK_ERR( fei->loadNodeBCs(numBCNodes, nodeIDs, fieldID,
00917           &offsets[0], values) );
00918 
00919   return(0);
00920 }
00921 
00922 //==============================================================================
00923 int init_elem_connectivities(fei::MatrixGraph* matrixGraph,
00924            PoissonData& poissonData)
00925 {
00926   //first load the information that defines this element block, and
00927   //the topology of each element in this element block.
00928 
00929   GlobalID elemBlockID = poissonData.getElemBlockID();
00930   int numLocalElements = poissonData.getNumLocalElements();
00931   int numNodesPerElement = poissonData.getNumNodesPerElement();
00932   int** fieldIDsTable = poissonData.getNodalFieldIDsTable();
00933 
00934   int nodeIDType = 0;
00935 
00936   int patternID =
00937     matrixGraph->definePattern(numNodesPerElement,
00938            nodeIDType, fieldIDsTable[0][0]);
00939 
00940   CHK_ERR( matrixGraph->initConnectivityBlock(elemBlockID,
00941               numLocalElements, patternID) );
00942 
00943   //now let's loop over all of the local elements, giving their 
00944   //nodal connectivity lists to the matrixGraph object.
00945 
00946   GlobalID* elemIDs = poissonData.getLocalElementIDs();
00947 
00948   for(int elem=0; elem<numLocalElements; elem++) {
00949     GlobalID* elemConnectivity =
00950       poissonData.getElementConnectivity(elemIDs[elem]);
00951 
00952     CHK_ERR( matrixGraph->initConnectivity(elemBlockID, elemIDs[elem],
00953            elemConnectivity) );
00954   }
00955 
00956   return(0);
00957 }
00958 
00959 //==============================================================================
00960 int set_shared_nodes(fei::VectorSpace* nodeSpace, PoissonData& poissonData)
00961 {
00962    int numLocalElements = poissonData.getNumLocalElements();
00963    int maxNumSharedNodes = (int)std::sqrt((double)numLocalElements);
00964    GlobalID* sharedNodeIDs = new GlobalID[maxNumSharedNodes];
00965    int* numProcsPerSharedNode = new int[maxNumSharedNodes];
00966    int** sharingProcs = new int*[maxNumSharedNodes];
00967    for(int i=0; i<maxNumSharedNodes; i++) sharingProcs[i] = new int[4];
00968 
00969    int numShared;
00970 
00971    //first, get the shared-node data for the left edge of the local block
00972 
00973    poissonData.getLeftSharedNodes(numShared, sharedNodeIDs,
00974                                   numProcsPerSharedNode, sharingProcs);
00975    int nodeIDType = 0;
00976 
00977    CHK_ERR( nodeSpace->initSharedIDs(numShared, nodeIDType, sharedNodeIDs,
00978                                 numProcsPerSharedNode, sharingProcs));
00979 
00980    //now, get the shared-node data for the right edge of the local block
00981 
00982    poissonData.getRightSharedNodes(numShared, sharedNodeIDs,
00983                                   numProcsPerSharedNode, sharingProcs);
00984 
00985    CHK_ERR( nodeSpace->initSharedIDs(numShared, nodeIDType, sharedNodeIDs,
00986                                 numProcsPerSharedNode, sharingProcs));
00987 
00988    //now, get the shared-node data for the bottom edge of the local block
00989 
00990    poissonData.getBottomSharedNodes(numShared, sharedNodeIDs,
00991                                   numProcsPerSharedNode, sharingProcs);
00992 
00993    CHK_ERR( nodeSpace->initSharedIDs(numShared, nodeIDType, sharedNodeIDs,
00994                                 numProcsPerSharedNode, sharingProcs));
00995 
00996    //finally, get the shared-node data for the top edge of the local block
00997 
00998    poissonData.getTopSharedNodes(numShared, sharedNodeIDs,
00999                                   numProcsPerSharedNode, sharingProcs);
01000 
01001    CHK_ERR( nodeSpace->initSharedIDs(numShared, nodeIDType, sharedNodeIDs,
01002                                 numProcsPerSharedNode, sharingProcs));
01003 
01004    for(int j=0; j<maxNumSharedNodes; j++) delete [] sharingProcs[j];
01005    delete [] sharingProcs;
01006    delete [] numProcsPerSharedNode;
01007    delete [] sharedNodeIDs;
01008 
01009    return(0);
01010 }
01011 
01012 //==============================================================================
01013 int load_elem_data(fei::MatrixGraph* matrixGraph,
01014        fei::Matrix* mat, fei::Vector* rhs,
01015        PoissonData& poissonData)
01016 {
01017   GlobalID elemBlockID = poissonData.getElemBlockID();
01018   int numLocalElements = poissonData.getNumLocalElements();
01019   GlobalID* elemIDs = poissonData.getLocalElementIDs();
01020 
01021   int numIndices = matrixGraph->getConnectivityNumIndices(elemBlockID);
01022 
01023   std::vector<int> indicesArray(numIndices);
01024   int* indicesPtr = &indicesArray[0];
01025 
01026   for(int elem=0; elem<numLocalElements; elem++) {
01027     double** elemStiffness = poissonData.getElemStiffness(elemIDs[elem]);
01028 
01029     int checkNumIndices = 0;
01030     CHK_ERR( matrixGraph->getConnectivityIndices(elemBlockID, elemIDs[elem],
01031                  numIndices, indicesPtr,
01032                  checkNumIndices) );
01033     if (checkNumIndices != numIndices) return(-1);
01034 
01035     CHK_ERR( mat->sumIn(elemBlockID, elemIDs[elem],
01036       elemStiffness));
01037 
01038     double* elemLoad = poissonData.getElemLoad(elemIDs[elem]);
01039 
01040     CHK_ERR( rhs->sumIn(numIndices, indicesPtr, elemLoad));
01041   }
01042 
01043   return(0);
01044 }
01045 
01046 //==============================================================================
01047 int load_BC_data(fei::LinearSystem* linSys, PoissonData& poissonData)
01048 {
01049   //first, have the data object generate the BC data
01050   poissonData.calculateBCs();
01051 
01052   int numBCNodes = poissonData.getNumBCNodes();
01053   GlobalID* nodeIDs = poissonData.getBCNodeIDs();
01054   int fieldID = poissonData.getBCFieldID();
01055   double* values = poissonData.getBCValues();
01056 
01057   CHK_ERR( linSys->loadEssentialBCs(numBCNodes, nodeIDs, 0, fieldID, 0, values) );
01058 
01059   return(0);
01060 }

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