FEI Version of the Day
Poisson_Elem.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 #include <fei_iostream.hpp>
00044 #include <fei_defs.h>
00045 #include <test_utils/Poisson_Elem.hpp>
00046 
00047 #include <cmath>
00048 #include <cstdlib>
00049 #include <assert.h>
00050 
00051 //==============================================================================
00052 Poisson_Elem::Poisson_Elem() {
00053 
00054     globalElemID_ = (GlobalID)(-1);
00055     ID_IsSet_ = false;
00056     numElemNodes_ = 4;
00057     numElemRows_ = 4;
00058 
00059     nodeList_ = NULL;
00060     nodalX_ = NULL;
00061     nodalY_ = NULL;
00062 
00063     elemStiff_ = NULL;
00064     elemLoad_ = NULL;
00065     internalsAllocated_ = false;
00066 
00067     elemLength_ = 0.0;
00068     elemLengthIsSet_ = false;
00069     totalLength_ = 0.0;
00070     totalLengthIsSet_ = false;
00071 
00072     loadAllocated_ = false;
00073     stiffAllocated_ = false;
00074 }
00075 
00076 //==============================================================================
00077 Poisson_Elem::~Poisson_Elem() {
00078    deleteMemory();
00079 }
00080 
00081 //==============================================================================
00082 void Poisson_Elem::deleteMemory()
00083 {
00084     if (!internalsAllocated_) return;
00085 
00086     delete [] nodeList_;
00087     delete [] nodalX_;
00088     delete [] nodalY_;
00089 
00090     internalsAllocated_ = false;
00091 
00092     if (loadAllocated_) {
00093        delete [] elemLoad_;
00094        loadAllocated_ = false;
00095     }
00096 
00097     if (stiffAllocated_) {
00098        delete [] elemStiff_[0];
00099        delete [] elemStiff_;
00100        stiffAllocated_ = false;
00101    }
00102 }
00103 
00104 //==============================================================================
00105 int Poisson_Elem::allocateInternals(int DOF) {
00106 
00107     assert(DOF==1);
00108 
00109     nodeList_ = new GlobalID[numElemNodes_];
00110     if (!nodeList_) return(1);
00111 
00112     nodalX_ = new double[numElemNodes_];
00113     if (!nodalX_) return(1);
00114 
00115     nodalY_ = new double[numElemNodes_];
00116     if (!nodalY_) return(1);
00117 
00118     for(int i = 0; i < numElemNodes_; i++) {
00119         nodeList_[i] = (GlobalID)0;
00120         nodalX_[i] = 0.0;
00121         nodalY_[i] = 0.0;
00122     }
00123 
00124     internalsAllocated_ = true;
00125     return(0);
00126 }
00127 
00128 //==============================================================================
00129 int Poisson_Elem::allocateLoad(int DOF)
00130 {
00131     assert(DOF==1);
00132 
00133     elemLoad_ = new double[numElemNodes_];
00134     if (!elemLoad_) return(1);
00135 
00136     loadAllocated_ = true;
00137     return(0);
00138 }
00139 
00140 //==============================================================================
00141 int Poisson_Elem::allocateStiffness(int DOF)
00142 {
00143     assert(DOF==1);
00144 
00145     elemStiff_ = new double*[numElemNodes_];
00146     if (!elemStiff_) return(1);
00147 
00148     elemStiff_[0] = NULL;
00149     elemStiff_[0] = new double[numElemNodes_*numElemNodes_];
00150     if (!elemStiff_[0]) return(1);
00151 
00152     int i;
00153     for(i=0; i<numElemNodes_*numElemNodes_; i++) elemStiff_[0][i] = 0.0;
00154 
00155     for(i=1; i<numElemNodes_; i++) {
00156         elemStiff_[i] = elemStiff_[i-1] + numElemNodes_;
00157     }
00158 
00159     stiffAllocated_ = true;
00160     return(0);
00161 }
00162 
00163 //==============================================================================
00164 GlobalID* Poisson_Elem::getElemConnPtr(int& size) {
00165 
00166     if (internalsAllocated_) {
00167         size = numElemNodes_;
00168         return(nodeList_);
00169     }
00170     else {
00171         size = 0;
00172         return(NULL);
00173     }
00174 }
00175 
00176 //==============================================================================
00177 double* Poisson_Elem::getElemLoad(int& size) {
00178         
00179     if (internalsAllocated_) {
00180         size = numElemNodes_;
00181         return(elemLoad_);
00182     }   
00183     else {
00184         size = 0;
00185         return(NULL);
00186     }
00187 }
00188 
00189 //==============================================================================
00190 double** Poisson_Elem::getElemStiff(int& size) {
00191         
00192     if (internalsAllocated_) {
00193         size = numElemNodes_*numElemNodes_;
00194         return(elemStiff_);
00195     }   
00196     else {
00197         size = 0;
00198         return(NULL);
00199     }
00200 }
00201 
00202 //==============================================================================
00203 void Poisson_Elem::calculateCoords() {
00204 //
00205 //This function calculates nodal x- and y-coordinates for this element.
00206 //NOTE: element IDs are assumed to be 1-based.
00207 //
00208     if (!internalsAllocated_)
00209         messageAbort("calculateCoords: internals not allocated.");
00210     if (!elemLengthIsSet_)
00211         messageAbort("calculateCoords: elemLength not set.");
00212     if (!totalLengthIsSet_)
00213         messageAbort("calculateCoords: totalLength not set.");
00214     if (!ID_IsSet_)
00215         messageAbort("calculateCoords: elemID not set.");
00216     if (std::abs(elemLength_) < 1.e-49)
00217         messageAbort("calculateCoords: elemLength == 0.");
00218 
00219     int lowLeft = 0;
00220     int lowRight = 1;
00221     int upperRight = 2;
00222     int upperLeft = 3;
00223 
00224     int elemsPerSide = (int)std::ceil(totalLength_/elemLength_);
00225 
00226     int elemX = (int)globalElemID_%elemsPerSide;
00227     if (elemX==0) elemX = elemsPerSide;
00228 
00229     int elemY = ((int)globalElemID_ - elemX)/elemsPerSide + 1;
00230 
00231     //elemX and elemY are 1-based coordinates of this element in
00232     //the global square of elements. The origin is position (1,1),
00233     //which is at the bottom left of the square.
00234 
00235     nodalX_[lowLeft] = (elemX-1)*elemLength_;
00236     nodalX_[upperLeft] = nodalX_[lowLeft];
00237     nodalX_[lowRight] = elemX*elemLength_;
00238     nodalX_[upperRight] = nodalX_[lowRight];
00239 
00240     nodalY_[lowLeft] = (elemY-1)*elemLength_;
00241     nodalY_[lowRight] = nodalY_[lowLeft];
00242     nodalY_[upperLeft] = elemY*elemLength_;
00243     nodalY_[upperRight] = nodalY_[upperLeft];
00244 }
00245 
00246 //==============================================================================
00247 void Poisson_Elem::messageAbort(const char* str) {
00248     fei::console_out() << "Poisson_Elem: ERROR: " << str << " Aborting." << FEI_ENDL;
00249     std::abort();
00250 }
00251 
00252 //==============================================================================
00253 void Poisson_Elem::calculateLoad() {
00254 
00255     assert (numElemRows_ == 4);
00256     
00257     elemLoad_[0] = -2.0*elemLength_*elemLength_;
00258     elemLoad_[1] = -2.0*elemLength_*elemLength_;
00259     elemLoad_[2] = -2.0*elemLength_*elemLength_;
00260     elemLoad_[3] = -2.0*elemLength_*elemLength_;
00261 }
00262 
00263 //==============================================================================
00264 void Poisson_Elem::calculateStiffness() {
00265 
00266     assert (numElemRows_ == 4);
00267     
00268     elemStiff_[0][0] = 2.0;
00269     elemStiff_[0][1] = -1.0;
00270     elemStiff_[0][2] = 0.0;
00271     elemStiff_[0][3] = -1.0;
00272 
00273     elemStiff_[1][0] = -1.0;
00274     elemStiff_[1][1] = 2.0;
00275     elemStiff_[1][2] = -1.0;
00276     elemStiff_[1][3] = 0.0;
00277 
00278     elemStiff_[2][0] = 0.0;
00279     elemStiff_[2][1] = -1.0;
00280     elemStiff_[2][2] = 2.0;
00281     elemStiff_[2][3] = -1.0;
00282 
00283     elemStiff_[3][0] = -1.0;
00284     elemStiff_[3][1] = 0.0;
00285     elemStiff_[3][2] = -1.0;
00286     elemStiff_[3][3] = 2.0;
00287 }
00288 
00289 
00290 
00291 //-----------------------------------------------------------------------
00292 
00293 void Poisson_Elem::dumpToScreen() {
00294 
00295     int i;
00296     
00297     FEI_COUT << " globalElemID_ = " << (int)globalElemID_ << FEI_ENDL;
00298     FEI_COUT << " numElemRows_  = " << numElemRows_ << FEI_ENDL;
00299     FEI_COUT << " elemLength_  = " << elemLength_ << FEI_ENDL;
00300     FEI_COUT << " the " << numElemNodes_ << " nodes: ";
00301     for (i = 0; i < numElemNodes_; ++i) {
00302       FEI_COUT << "   " << (int)nodeList_[i];
00303     }
00304     FEI_COUT << FEI_ENDL;
00305     FEI_COUT << " the " << numElemRows_ << " load vector terms: ";
00306     for (i = 0; i < numElemRows_; ++i) {
00307       FEI_COUT << "   " << elemLoad_[i];
00308     }
00309     FEI_COUT << FEI_ENDL << FEI_ENDL;
00310     return;
00311 }
00312 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends