Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/src/Discretization/Integration/Intrepid_CubatureTensorDef.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
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 Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00049 namespace Intrepid {
00050 
00051 template <class Scalar, class ArrayPoint, class ArrayWeight>
00052 CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::CubatureTensor(std::vector< Teuchos::RCP<Cubature<Scalar,ArrayPoint,ArrayWeight> > > cubatures) {
00053   unsigned numCubs = cubatures.size();
00054   TEUCHOS_TEST_FOR_EXCEPTION( (numCubs < 1),
00055                       std::out_of_range,
00056                       ">>> ERROR (CubatureTensor): Input cubature array must be of size 1 or larger.");
00057 
00058   cubatures_ = cubatures;
00059 
00060   std::vector<int> tmp;
00061   unsigned numDegrees = 0;
00062   for (unsigned i=0; i<numCubs; i++) {
00063     cubatures[i]->getAccuracy(tmp);
00064     numDegrees += tmp.size();
00065   }
00066 
00067   degree_.assign(numDegrees, 0);
00068   int count  = 0;
00069   dimension_ = 0;
00070   for (unsigned i=0; i<numCubs; i++) {
00071     cubatures[i]->getAccuracy(tmp);
00072     for (unsigned j=0; j<tmp.size(); j++) {
00073       degree_[count] = tmp[j];
00074       count++;
00075     }
00076     dimension_ += cubatures[i]->getDimension();
00077   }
00078 }
00079 
00080 
00081 
00082 template <class Scalar, class ArrayPoint, class ArrayWeight>
00083 CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::CubatureTensor(Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature1,
00084                                                               Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature2) {
00085   cubatures_.resize(2);
00086   cubatures_[0] = cubature1;
00087   cubatures_[1] = cubature2;
00088 
00089   degree_.assign(2, 0);
00090   std::vector<int> d(1);
00091   cubatures_[0]->getAccuracy(d); degree_[0] = d[0];
00092   cubatures_[1]->getAccuracy(d); degree_[1] = d[0];
00093 
00094   dimension_ = cubatures_[0]->getDimension() + cubatures_[1]->getDimension();
00095 }
00096 
00097 
00098 
00099 template <class Scalar, class ArrayPoint, class ArrayWeight>
00100 CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::CubatureTensor(Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature1,
00101                                                               Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature2,
00102                                                               Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature3) {
00103   cubatures_.resize(3);
00104   cubatures_[0] = cubature1;
00105   cubatures_[1] = cubature2;
00106   cubatures_[2] = cubature3;
00107 
00108   degree_.assign(3, 0);
00109   std::vector<int> d(1);
00110   cubatures_[0]->getAccuracy(d); degree_[0] = d[0];
00111   cubatures_[1]->getAccuracy(d); degree_[1] = d[0];
00112   cubatures_[2]->getAccuracy(d); degree_[2] = d[0];
00113 
00114   dimension_ = cubatures_[0]->getDimension() + cubatures_[1]->getDimension() + cubatures_[2]->getDimension();
00115 }
00116 
00117 
00118 
00119 template <class Scalar, class ArrayPoint, class ArrayWeight>
00120 CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::CubatureTensor(Teuchos::RCP<CubatureDirect<Scalar,ArrayPoint,ArrayWeight> > cubature, int n) {
00121   cubatures_.resize(n);
00122   for (int i=0; i<n; i++) {
00123     cubatures_[i] = cubature;
00124   }
00125 
00126   std::vector<int> d(1);
00127   cubatures_[0]->getAccuracy(d);
00128   degree_.assign(n,d[0]);
00129 
00130   dimension_ = cubatures_[0]->getDimension()*n;
00131 }
00132 
00133 
00134 
00135 template <class Scalar, class ArrayPoint, class ArrayWeight>
00136 void CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::getCubature(ArrayPoint  & cubPoints,
00137                                                                 ArrayWeight & cubWeights) const {
00138   int numCubPoints = getNumPoints();
00139   int cubDim       = getDimension();
00140   // check size of cubPoints and cubWeights
00141   TEUCHOS_TEST_FOR_EXCEPTION( ( ( (int)cubPoints.size() < numCubPoints*cubDim ) || ( (int)cubWeights.size() < numCubPoints ) ),
00142                       std::out_of_range,
00143                       ">>> ERROR (CubatureTensor): Insufficient space allocated for cubature points or weights.");
00144 
00145   unsigned numCubs   = cubatures_.size();
00146   std::vector<unsigned> numLocPoints(numCubs);
00147   std::vector<unsigned> locDim(numCubs);
00148   std::vector< FieldContainer<Scalar> > points(numCubs);
00149   std::vector< FieldContainer<Scalar> > weights(numCubs);
00150 
00151   // extract required points and weights
00152   for (unsigned i=0; i<numCubs; i++) {
00153 
00154     numLocPoints[i] = cubatures_[i]->getNumPoints();
00155     locDim[i]       = cubatures_[i]->getDimension();
00156     points[i].resize(numLocPoints[i], locDim[i]);
00157     weights[i].resize(numLocPoints[i]);
00158 
00159     // cubPoints and cubWeights are used here only for temporary data retrieval
00160     cubatures_[i]->getCubature(cubPoints, cubWeights);
00161     for (unsigned pt=0; pt<numLocPoints[i]; pt++) {
00162       for (unsigned d=0; d<locDim[i]; d++) {
00163         points[i](pt,d) = cubPoints(pt,d);
00164         weights[i](pt)  = cubWeights(pt);
00165       }
00166     }
00167 
00168   }
00169 
00170   // reset all weights to 1.0
00171   for (int i=0; i<numCubPoints; i++) {
00172       cubWeights(i) = (Scalar)1.0;
00173   }
00174 
00175   // fill tensor-product cubature
00176   int globDimCounter = 0;
00177   int shift          = 1;
00178   for (unsigned i=0; i<numCubs; i++) {
00179 
00180     for (int j=0; j<numCubPoints; j++) {
00181       /* int itmp = ((j*shift) % numCubPoints) + (j / (numCubPoints/shift)); // equivalent, but numerically unstable */
00182       int itmp = (j % (numCubPoints/shift))*shift + (j / (numCubPoints/shift));
00183       for (unsigned k=0; k<locDim[i]; k++) {
00184         cubPoints(itmp , globDimCounter+k) = points[i](j % numLocPoints[i], k);
00185       }
00186       cubWeights( itmp ) *= weights[i](j % numLocPoints[i]);
00187     }
00188     
00189     shift *= numLocPoints[i];
00190     globDimCounter += locDim[i];
00191   }
00192 
00193 } // end getCubature
00194 
00195 
00196 
00197 template <class Scalar, class ArrayPoint, class ArrayWeight>
00198 int CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::getNumPoints() const {
00199   unsigned numCubs = cubatures_.size();
00200   int numCubPoints = 1;
00201   for (unsigned i=0; i<numCubs; i++) {
00202     numCubPoints *= cubatures_[i]->getNumPoints();
00203   }
00204   return numCubPoints;
00205 } // end getNumPoints
00206 
00207 
00208 template <class Scalar, class ArrayPoint, class ArrayWeight>
00209 int CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::getDimension() const {
00210   return dimension_;
00211 } // end dimension
00212 
00213 
00214 
00215 template <class Scalar, class ArrayPoint, class ArrayWeight>
00216 void CubatureTensor<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(std::vector<int> & degree) const {
00217   degree = degree_;
00218 } // end getAccuracy
00219 
00220 } // end namespace Intrepid