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