00001 /* @HEADER@ */ 00002 // ************************************************************************ 00003 // 00004 // Sundance 00005 // Copyright (2005) Sandia Corporation 00006 // 00007 // Copyright (year first published) Sandia Corporation. Under the terms 00008 // of Contract DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 00009 // retains certain rights in this software. 00010 // 00011 // This library is free software; you can redistribute it and/or modify 00012 // it under the terms of the GNU Lesser General Public License as 00013 // published by the Free Software Foundation; either version 2.1 of the 00014 // License, or (at your option) any later version. 00015 // 00016 // This library is distributed in the hope that it will be useful, but 00017 // WITHOUT ANY WARRANTY; without even the implied warranty of 00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00019 // Lesser General Public License for more details. 00020 // 00021 // You should have received a copy of the GNU Lesser General Public 00022 // License along with this library; if not, write to the Free Software 00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00024 // USA 00025 // Questions? Contact Kevin Long (krlong@sandia.gov), 00026 // Sandia National Laboratories, Livermore, California, USA 00027 // 00028 // ************************************************************************ 00029 /* @HEADER@ */ 00030 00031 #ifndef SUNDANCE_CELLJACOBIANBATCH_H 00032 #define SUNDANCE_CELLJACOBIANBATCH_H 00033 00034 00035 #include "SundanceDefs.hpp" 00036 #include "SundanceObjectWithVerbosity.hpp" 00037 #include "Teuchos_Array.hpp" 00038 00039 namespace Sundance 00040 { 00041 using namespace Teuchos; 00042 00043 00044 00045 00046 /** 00047 * A CellJacobianBatch is a collection of Jacobian matrices 00048 * for many quadrature 00049 * points distributed over batch of cells. All cells 00050 * must have the same dimension 00051 * and number of quadrature points. Affine cells have constant Jacobians, 00052 * so in that case the quadrature point index can be ignored. 00053 * See the ReferenceCellBase documentation 00054 * for definitions of the coordinate systems used. 00055 * 00056 * <H4> Data layout </H4> 00057 * All Jacobian elements for all points and cells are 00058 * packed into a single vector. 00059 * The length of the vector is 00060 * \f$ N_{cell} \times N_{quad} \times D^2, \f$ where 00061 * \f$D\f$ is the spatial dimension. 00062 * The indices into the vector cycle in the following order 00063 * (from slowest to fastest) 00064 * <ol> 00065 * <li> cell number 00066 * <li> quadrature point number 00067 * <li> physical coordinate direction 00068 * <li> reference coordinate direction 00069 * </ol> 00070 * Thus, the jacobian values for the \f$q\f$-th quadrature point on the 00071 * \f$c\f$-th cell are the \f$D^2\f$ entries following 00072 * the \f$(q + cN_{quad})D^2\f$-th 00073 * element. 00074 * 00075 */ 00076 class CellJacobianBatch 00077 : public ObjectWithClassVerbosity<CellJacobianBatch> 00078 { 00079 public: 00080 /** empty ctor */ 00081 CellJacobianBatch(); 00082 00083 /** get the spatial dimension */ 00084 int spatialDim() const {return spatialDim_;} 00085 00086 /** get the cell dimension */ 00087 int cellDim() const {return cellDim_;} 00088 00089 /** get the number of cells in the batch */ 00090 int numCells() const {return numCells_;} 00091 00092 /** get the number of quad points per cell */ 00093 int numQuadPoints() const {return numQuad_;} 00094 00095 /** resize the batch */ 00096 void resize(int numCells, int numQuad, int spatialDim, int cellDim); 00097 00098 /** resize the batch, using one quadrature point per cell 00099 * (appropriate for affine elements) */ 00100 void resize(int numCells, int spatialDim, int cellDim); 00101 00102 /** Get a pointer to the values at the q-th quadrature 00103 * point on the c-th cell. 00104 * @param c the index of the cell in the batch 00105 * @param q the index of the quadrature point 00106 */ 00107 double* jVals(int c, int q); 00108 00109 /** 00110 * Get a pointer to the start of the c-th Jacobian in the batch. 00111 */ 00112 double* jVals(int c) 00113 {return &(J_[c*jSize_]);} 00114 00115 /** Get a constant pointer to start of c-th Jacobian in the batch */ 00116 const double *jVals(int c) const { return &(J_[c*jSize_]); } 00117 00118 /** */ 00119 double* detJ(int c) 00120 {return &(detJ_[c]);} 00121 00122 /** get the vector of determinant values */ 00123 const Array<double>& detJ() const 00124 {if (!isFactored_ && cellDim()==spatialDim()) factor(); return detJ_;} 00125 00126 /** 00127 * Apply a cell's inverse Jacobian to (possibly) multiple rhs 00128 * stored in column-major order. 00129 */ 00130 void applyInvJ(int cell, int q, double* rhs, 00131 int nRhs, bool trans) const ; 00132 00133 /** 00134 * Apply an affine cell's inverse Jacobian to (possibly) multiple rhs 00135 * stored in column-major order. 00136 */ 00137 void applyInvJ(int cell, double* rhs, 00138 int nRhs, bool trans) const 00139 {applyInvJ(cell, 0, rhs, nRhs, trans);} 00140 00141 /** 00142 * Get the explicit inverse of the Jacobian for the given 00143 * (cell, quad) combination. 00144 */ 00145 void getInvJ(int cell, int quad, Array<double>& invJ) const ; 00146 00147 /** 00148 * Get the explicit inverse of the Jacobian for the given 00149 * affine cell. 00150 */ 00151 void getInvJ(int cell, Array<double>& invJ) const 00152 {getInvJ(cell, 0, invJ);} 00153 00154 00155 /** */ 00156 void print(std::ostream& os) const ; 00157 00158 static double& totalFlops() {static double rtn = 0; return rtn;} 00159 00160 00161 00162 static void addFlops(const double& flops) {totalFlops() += flops;} 00163 00164 private: 00165 00166 void factor() const ; 00167 00168 void computeInverses() const ; 00169 00170 int spatialDim_; 00171 int cellDim_; 00172 int jSize_; 00173 int numCells_; 00174 int numQuad_; 00175 mutable Array<int> iPiv_; 00176 mutable Array<double> J_; 00177 mutable Array<double> detJ_; 00178 mutable Array<double> invJ_; 00179 mutable bool isFactored_; 00180 mutable bool hasInverses_; 00181 mutable bool hasDetJ_; 00182 }; 00183 00184 00185 inline std::ostream& operator<<(std::ostream& os, 00186 const CellJacobianBatch& J) 00187 { 00188 J.print(os); 00189 return os; 00190 } 00191 } 00192 00193 00194 00195 #endif