00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceCellJacobianBatch.hpp"
00032 #include "PlayaExceptions.hpp"
00033 #include "SundanceOut.hpp"
00034 #include "PlayaTabs.hpp"
00035 #include "Teuchos_Time.hpp"
00036 #include "Teuchos_TimeMonitor.hpp"
00037
00038 using namespace Sundance;
00039 using namespace Sundance;
00040 using namespace Teuchos;
00041 using namespace Sundance;
00042
00043
00044
00045 extern "C"
00046 {
00047
00048 void dgetrs_(const char* trans, const int* N, const int* NRHS,
00049 const double* A, const int* lda,
00050 const int* iPiv, double* B, const int* ldb, int* info);
00051
00052
00053 void dgetrf_(const int* M, const int* N, double* A, const int* lda,
00054 const int* iPiv, int* info);
00055 }
00056
00057 static Time& jacobianInversionTimer()
00058 {
00059 static RCP<Time> rtn
00060 = TimeMonitor::getNewTimer("jacobian inversion");
00061 return *rtn;
00062 }
00063
00064 static Time& jacobianFactoringTimer()
00065 {
00066 static RCP<Time> rtn
00067 = TimeMonitor::getNewTimer("jacobian factoring");
00068 return *rtn;
00069 }
00070
00071
00072 CellJacobianBatch::CellJacobianBatch()
00073 : spatialDim_(0), cellDim_(0),
00074 jSize_(0), numCells_(0), numQuad_(0), iPiv_(), J_(), detJ_(), invJ_(),
00075 isFactored_(false), hasInverses_(false)
00076 {}
00077
00078 void CellJacobianBatch::resize(int numCells, int numQuad,
00079 int spatialDim, int cellDim)
00080 {
00081 spatialDim_ = spatialDim;
00082 cellDim_ = cellDim;
00083 if (spatialDim_ == cellDim_)
00084 {
00085 jSize_ = spatialDim_*spatialDim_;
00086 }
00087
00088 numCells_ = numCells;
00089 numQuad_ = numQuad;
00090 iPiv_.resize(spatialDim_*numCells_*numQuad_);
00091 J_.resize(spatialDim_*spatialDim_*numCells_*numQuad_);
00092 detJ_.resize(numCells_*numQuad_);
00093 isFactored_ = false;
00094 hasInverses_ = false;
00095 }
00096
00097 void CellJacobianBatch::resize(int numCells, int spatialDim, int cellDim)
00098 {
00099 spatialDim_ = spatialDim;
00100 cellDim_ = cellDim;
00101 if (spatialDim_ == cellDim_)
00102 {
00103 jSize_ = spatialDim_*spatialDim_;
00104 }
00105
00106 numCells_ = numCells;
00107 numQuad_ = 1;
00108 iPiv_.resize(spatialDim_*numCells_);
00109 J_.resize(spatialDim_*spatialDim_*numCells_);
00110 detJ_.resize(numCells_);
00111 isFactored_ = false;
00112 hasInverses_ = false;
00113 }
00114
00115 void CellJacobianBatch::factor() const
00116 {
00117 TimeMonitor timer(jacobianFactoringTimer());
00118 if (isFactored_) return;
00119
00120
00121
00122
00123
00124
00125 TEUCHOS_TEST_FOR_EXCEPTION(spatialDim_ != cellDim_, std::logic_error,
00126 "Attempting to factor the Jacobian of a cell "
00127 "that is not of maximal dimension");
00128 Tabs tabs;
00129 SUNDANCE_OUT(this->verb() > 2,
00130 tabs << "factoring Jacobians");
00131
00132 for (int cell=0; cell<numCells_; cell++)
00133 {
00134 for (int q=0; q<numQuad_; q++)
00135 {
00136 int start = (cell*numQuad_ + q)*jSize_;
00137
00138
00139 double* jFactPtr = &(J_[start]);
00140 int* iPiv = &(iPiv_[(q + cell*numQuad_)*spatialDim_]);
00141
00142
00143 int lda = spatialDim_;
00144
00145 int info = 0;
00146
00147
00148 ::dgetrf_( &spatialDim_, &spatialDim_, jFactPtr, &lda, iPiv, &info);
00149
00150 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
00151 "CellJacobianBatch::setJacobian(): factoring failed");
00152
00153
00154
00155 double detJ = 1.0;
00156 for (int i=0; i<spatialDim_; i++)
00157 {
00158 detJ *= jFactPtr[i + spatialDim_*i];
00159 }
00160 detJ_[cell*numQuad_ + q] = detJ;
00161 }
00162 }
00163
00164 addFlops(numCells_ * spatialDim_ * (1.0 + spatialDim_ * spatialDim_));
00165 isFactored_ = true;
00166 }
00167
00168 void CellJacobianBatch::computeInverses() const
00169 {
00170 TimeMonitor timer(jacobianInversionTimer());
00171 if (hasInverses_) return;
00172
00173 Tabs tabs;
00174 SUNDANCE_OUT(this->verb() > 2,
00175 tabs << "inverting Jacobians");
00176
00177 invJ_.resize(spatialDim_*spatialDim_*numQuad_*numCells_);
00178
00179 if (!isFactored_) factor();
00180
00181 for (int cell=0; cell<numCells_; cell++)
00182 {
00183 for (int q=0; q<numQuad_; q++)
00184 {
00185 int start = (cell*numQuad_ + q)*jSize_;
00186
00187
00188 double* jFactPtr = &(J_[start]);
00189 double* invJPtr = &(invJ_[start]);
00190 int* iPiv = &(iPiv_[(q + cell*numQuad_)*spatialDim_]);
00191
00192 int info = 0;
00193
00194
00195 for (int i=0; i<spatialDim_; i++)
00196 {
00197 for (int j=0; j<spatialDim_; j++)
00198 {
00199 if (i==j) invJPtr[i*spatialDim_+j] = 1.0;
00200 else invJPtr[i*spatialDim_+j] = 0.0;
00201 }
00202 }
00203
00204 ::dgetrs_("N", &spatialDim_, &spatialDim_, jFactPtr,
00205 &spatialDim_, iPiv, invJPtr, &spatialDim_, &info);
00206
00207 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
00208 "CellJacobianBatch::setJacobian(): inversion failed");
00209 }
00210 }
00211 addFlops(numCells_ * spatialDim_ * spatialDim_);
00212 hasInverses_ = true;
00213 }
00214
00215 void CellJacobianBatch::applyInvJ(int cell, int q,
00216 double* rhs, int nRhs, bool trans) const
00217 {
00218 if (!isFactored_) factor();
00219
00220 double* jFactPtr = &(J_[(cell*numQuad_ + q)*spatialDim_*spatialDim_]);
00221 int* iPiv = &(iPiv_[(q + cell*numQuad_)*spatialDim_]);
00222
00223 int info = 0;
00224
00225 if (trans)
00226 {
00227 ::dgetrs_("T", &spatialDim_, &nRhs, jFactPtr, &spatialDim_,
00228 iPiv, rhs, &spatialDim_, &info);
00229 }
00230 else
00231 {
00232 ::dgetrs_("N", &spatialDim_, &nRhs, jFactPtr, &spatialDim_,
00233 iPiv, rhs, &spatialDim_, &info);
00234 }
00235
00236 addFlops(numCells_ * spatialDim_ * spatialDim_ * nRhs);
00237 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::runtime_error,
00238 "CellJacobianBatch::applyInvJ(): backsolve failed");
00239 }
00240
00241 void CellJacobianBatch::getInvJ(int cell, int quad, Array<double>& invJ) const
00242 {
00243 if (!hasInverses_) computeInverses();
00244
00245 int start = (cell*numQuad_ + quad)*jSize_;
00246
00247 invJ.resize(spatialDim_*spatialDim_);
00248
00249 for (int col=0; col<spatialDim_; col++)
00250 {
00251 for (int row=0; row<spatialDim_; row++)
00252 {
00253 invJ[col + spatialDim_*row] = invJ_[start + col + spatialDim_*row];
00254 }
00255 }
00256 }
00257
00258 void CellJacobianBatch::print(std::ostream& os) const
00259 {
00260 if (!hasInverses_) computeInverses();
00261
00262 for (int c=0; c<numCells_; c++)
00263 {
00264 os << "cell " << c << std::endl;
00265 for (int q=0; q<numQuad_; q++)
00266 {
00267 int start = (c*numQuad_ + q)*jSize_;
00268 if (numQuad_ > 1) os << "q=" << q << " ";
00269 os << "{";
00270 for (int i=0; i<spatialDim_; i++)
00271 {
00272 if (i != 0) os << ", ";
00273 os << "{";
00274 for (int j=0; j<spatialDim_; j++)
00275 {
00276 if (j != 0) os << ", ";
00277 os << invJ_[start + i*spatialDim_ + j];
00278 }
00279 os << "}";
00280 }
00281 os << "}" << std::endl;
00282 }
00283
00284 }
00285 }
00286
00287