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 "SundanceIntegralGroup.hpp"
00032 #include "SundanceRefIntegral.hpp"
00033 #include "SundanceQuadratureIntegral.hpp"
00034 #include "SundanceMaximalQuadratureIntegral.hpp"
00035 #include "SundanceCurveQuadratureIntegral.hpp"
00036 #include "SundanceOut.hpp"
00037 #include "PlayaTabs.hpp"
00038 #include "Teuchos_TimeMonitor.hpp"
00039
00040 using namespace Sundance;
00041 using namespace Teuchos;
00042
00043
00044 static Time& integrationTimer()
00045 {
00046 static RCP<Time> rtn
00047 = TimeMonitor::getNewTimer("integration");
00048 return *rtn;
00049 }
00050
00051
00052 IntegralGroup
00053 ::IntegralGroup(const Array<RCP<ElementIntegral> >& integrals,
00054 const Array<int>& resultIndices,
00055 int verb)
00056 : integrationVerb_(findIntegrationVerb(integrals)),
00057 transformVerb_(findTransformVerb(integrals)),
00058 order_(0),
00059 nTestNodes_(0),
00060 nUnkNodes_(0),
00061 testID_(),
00062 unkID_(),
00063 testBlock_(),
00064 unkBlock_(),
00065 mvIndices_(),
00066 integrals_(integrals),
00067 resultIndices_(resultIndices),
00068 termUsesMaximalCofacets_(integrals_.size()),
00069 requiresMaximalCofacet_(SomeTermsNeedCofacets),
00070 derivs_()
00071 {
00072 Tabs tab;
00073 SUNDANCE_MSG2(verb, tab << "forming 0-form integral group");
00074 bool allReqMaximalCofacets = true;
00075 bool noneReqMaximalCofacets = true;
00076
00077
00078 for (int i=0; i<integrals_.size(); i++)
00079 {
00080 Tabs tab1;
00081 SUNDANCE_MSG2(verb, tab1 << "integral #" << i << ", numFacetCases="
00082 << integrals[i]->nFacetCases());
00083 if (integrals[i]->nFacetCases() > 1)
00084 {
00085 Tabs tab2;
00086
00087 noneReqMaximalCofacets = false;
00088 termUsesMaximalCofacets_[i] = true;
00089 SUNDANCE_MSG2(verb, tab2 << "I need maximal cofacets");
00090 }
00091 else
00092 {
00093 Tabs tab2;
00094 allReqMaximalCofacets = false;
00095 termUsesMaximalCofacets_[i] = false;
00096 SUNDANCE_MSG2(verb, tab2 << "I do not need maximal cofacets");
00097 }
00098 }
00099
00100 if (allReqMaximalCofacets)
00101 {
00102 requiresMaximalCofacet_ = AllTermsNeedCofacets;
00103 }
00104 else if (noneReqMaximalCofacets)
00105 {
00106 requiresMaximalCofacet_ = NoTermsNeedCofacets;
00107 }
00108
00109 Tabs tab1;
00110 SUNDANCE_MSG2(verb, tab1 << "result=" << requiresMaximalCofacet_);
00111 }
00112
00113
00114
00115
00116 IntegralGroup
00117 ::IntegralGroup(const Array<int>& testID,
00118 const Array<int>& testBlock,
00119 const Array<int>& mvIndices,
00120 const Array<RCP<ElementIntegral> >& integrals,
00121 const Array<int>& resultIndices,
00122 const Array<MultipleDeriv>& derivs,
00123 int verb)
00124 : integrationVerb_(findIntegrationVerb(integrals)),
00125 transformVerb_(findTransformVerb(integrals)),
00126 order_(1),
00127 nTestNodes_(0),
00128 nUnkNodes_(0),
00129 testID_(testID),
00130 unkID_(),
00131 testBlock_(testBlock),
00132 unkBlock_(),
00133 mvIndices_(mvIndices),
00134 integrals_(integrals),
00135 resultIndices_(resultIndices),
00136 termUsesMaximalCofacets_(integrals_.size()),
00137 requiresMaximalCofacet_(SomeTermsNeedCofacets),
00138 derivs_(derivs)
00139 {
00140 Tabs tab;
00141 SUNDANCE_MSG2(verb, tab << "forming 1-form integral group");
00142 bool allReqMaximalCofacets = true;
00143 bool noneReqMaximalCofacets = true;
00144
00145
00146 for (int i=0; i<integrals.size(); i++)
00147 {
00148 int nt = integrals[i]->nNodesTest();
00149 if (i > 0)
00150 {
00151 TEUCHOS_TEST_FOR_EXCEPTION(nt != nTestNodes_, std::logic_error,
00152 "IntegralGroup ctor detected integrals with inconsistent numbers of test nodes");
00153 }
00154 Tabs tab1;
00155 SUNDANCE_MSG2(verb, tab1 << "integral #" << i << ", numFacetCases="
00156 << integrals[i]->nFacetCases());
00157 nTestNodes_ = nt;
00158 if (integrals[i]->nFacetCases() > 1)
00159 {
00160 Tabs tab2;
00161 SUNDANCE_MSG2(verb, tab2 << "I need maximal cofacets");
00162
00163 noneReqMaximalCofacets = false;
00164 termUsesMaximalCofacets_[i] = true;
00165 }
00166 else
00167 {
00168 Tabs tab2;
00169 SUNDANCE_MSG2(verb, tab2 << "I do not need maximal cofacets");
00170 allReqMaximalCofacets = false;
00171 termUsesMaximalCofacets_[i] = false;
00172 }
00173 }
00174
00175 if (allReqMaximalCofacets)
00176 {
00177 requiresMaximalCofacet_ = AllTermsNeedCofacets;
00178 }
00179 else if (noneReqMaximalCofacets)
00180 {
00181 requiresMaximalCofacet_ = NoTermsNeedCofacets;
00182 }
00183
00184 Tabs tab1;
00185 SUNDANCE_MSG2(verb, tab1 << "result=" << requiresMaximalCofacet_);
00186 }
00187
00188 IntegralGroup
00189 ::IntegralGroup(const Array<int>& testID,
00190 const Array<int>& testBlock,
00191 const Array<int>& unkID,
00192 const Array<int>& unkBlock,
00193 const Array<RCP<ElementIntegral> >& integrals,
00194 const Array<int>& resultIndices,
00195 const Array<MultipleDeriv>& derivs,
00196 int verb)
00197 : integrationVerb_(findIntegrationVerb(integrals)),
00198 transformVerb_(findTransformVerb(integrals)),
00199 order_(2),
00200 nTestNodes_(0),
00201 nUnkNodes_(0),
00202 testID_(testID),
00203 unkID_(unkID),
00204 testBlock_(testBlock),
00205 unkBlock_(unkBlock),
00206 mvIndices_(),
00207 integrals_(integrals),
00208 resultIndices_(resultIndices),
00209 termUsesMaximalCofacets_(integrals_.size()),
00210 requiresMaximalCofacet_(SomeTermsNeedCofacets),
00211 derivs_(derivs)
00212 {
00213 Tabs tab;
00214 SUNDANCE_MSG2(verb, tab << "forming 2-form integral group");
00215 bool allReqMaximalCofacets = true;
00216 bool noneReqMaximalCofacets = true;
00217
00218
00219 for (int i=0; i<integrals.size(); i++)
00220 {
00221 Tabs tab1;
00222 SUNDANCE_MSG2(verb, tab1 << "integral #" << i << ", numFacetCases="
00223 << integrals[i]->nFacetCases());
00224 int nt = integrals[i]->nNodesTest();
00225 int nu = integrals[i]->nNodesUnk();
00226 if (i > 0)
00227 {
00228 TEUCHOS_TEST_FOR_EXCEPTION(nt != nTestNodes_, std::logic_error,
00229 "IntegralGroup ctor detected integrals with inconsistent numbers of test nodes");
00230 TEUCHOS_TEST_FOR_EXCEPTION(nu != nUnkNodes_, std::logic_error,
00231 "IntegralGroup ctor detected integrals with inconsistent numbers of unk nodes");
00232 }
00233 nTestNodes_ = nt;
00234 nUnkNodes_ = nu;
00235
00236 if (integrals[i]->nFacetCases() > 1)
00237 {
00238
00239 noneReqMaximalCofacets = false;
00240 termUsesMaximalCofacets_[i] = true;
00241 Tabs tab2;
00242 SUNDANCE_MSG2(verb, tab2 << "I need maximal cofacets");
00243 }
00244 else
00245 {
00246 allReqMaximalCofacets = false;
00247 termUsesMaximalCofacets_[i] = false;
00248 Tabs tab2;
00249 SUNDANCE_MSG2(verb, tab2 << "I do not need maximal cofacets");
00250 }
00251 }
00252
00253 if (allReqMaximalCofacets)
00254 {
00255 requiresMaximalCofacet_ = AllTermsNeedCofacets;
00256 }
00257 else if (noneReqMaximalCofacets)
00258 {
00259 requiresMaximalCofacet_ = NoTermsNeedCofacets;
00260 }
00261
00262 Tabs tab1;
00263 SUNDANCE_MSG2(verb, tab1 << "result=" << requiresMaximalCofacet_);
00264 }
00265
00266
00267 bool IntegralGroup
00268 ::evaluate(const CellJacobianBatch& JTrans,
00269 const CellJacobianBatch& JVol,
00270 const Array<int>& isLocalFlag,
00271 const Array<int>& facetIndex,
00272 const RCP<Array<int> >& cellLIDs,
00273 const Array<RCP<EvalVector> >& vectorCoeffs,
00274 const Array<double>& constantCoeffs,
00275 RCP<Array<double> >& A) const
00276 {
00277 TimeMonitor timer(integrationTimer());
00278 Tabs tab0(0);
00279
00280
00281 SUNDANCE_MSG1(integrationVerb(), tab0 << "evaluating integral group with "
00282 << integrals_.size() << " integrals");
00283
00284 SUNDANCE_MSG3(integrationVerb(),
00285 tab0 << "num integration cells = " << JVol.numCells());
00286 SUNDANCE_MSG3(integrationVerb(),
00287 tab0 << "num nodes in output = " << integrals_[0]->nNodes());
00288
00289
00290 if (integrals_[0]->nNodes() == -1) A->resize(1);
00291 else A->resize(JVol.numCells() * integrals_[0]->nNodes());
00292 double* aPtr = &((*A)[0]);
00293 int n = A->size();
00294 for (int i=0; i<n; i++) aPtr[i] = 0.0;
00295
00296 SUNDANCE_MSG5(integrationVerb(), tab0 << "begin A=");
00297 if (integrationVerb() >=5) writeTable(Out::os(), tab0, *A, 6);
00298
00299
00300 for (int i=0; i<integrals_.size(); i++)
00301 {
00302 Tabs tab1;
00303 SUNDANCE_MSG1(integrationVerb(), tab1 << "group member i=" << i
00304 << " of " << integrals_.size());
00305 Tabs tab2;
00306
00307 const RefIntegral* ref
00308 = dynamic_cast<const RefIntegral*>(integrals_[i].get());
00309 const QuadratureIntegral* quad
00310 = dynamic_cast<const QuadratureIntegral*>(integrals_[i].get());
00311 const MaximalQuadratureIntegral* maxQuad
00312 = dynamic_cast<const MaximalQuadratureIntegral*>(integrals_[i].get());
00313 const CurveQuadratureIntegral* curveQuad
00314 = dynamic_cast<const CurveQuadratureIntegral*>(integrals_[i].get());
00315
00316 if (ref!=0)
00317 {
00318 SUNDANCE_MSG2(integrationVerb(),
00319 tab2 << "Integrating term group " << i
00320 << " by transformed reference integral");
00321 double f = constantCoeffs[resultIndices_[i]];
00322 SUNDANCE_MSG2(integrationVerb(),
00323 tab2 << "Coefficient is " << f);
00324 ref->transform(JTrans, JVol, isLocalFlag, facetIndex, cellLIDs , f, A);
00325 }
00326 else if (quad != 0)
00327 {
00328 SUNDANCE_MSG2(integrationVerb(),
00329 tab2 << "Integrating term group " << i
00330 << " by quadrature");
00331
00332 TEUCHOS_TEST_FOR_EXCEPTION(vectorCoeffs[resultIndices_[i]]->length()==0,
00333 std::logic_error,
00334 "zero-length coeff vector detected in "
00335 "quadrature integration branch of "
00336 "IntegralGroup::evaluate(). std::string value is ["
00337 << vectorCoeffs[resultIndices_[i]]->str()
00338 << "]");
00339
00340 Tabs tab3;
00341 SUNDANCE_MSG3(integrationVerb(),
00342 tab3 << "coefficients are " << vectorCoeffs[resultIndices_[i]]->str());
00343
00344 const double* const f = vectorCoeffs[resultIndices_[i]]->start();
00345 quad->transform(JTrans, JVol, isLocalFlag, facetIndex, cellLIDs , f, A);
00346 }
00347 else if (maxQuad != 0)
00348 {
00349 SUNDANCE_MSG2(integrationVerb(),
00350 tab2 << "Integrating term group " << i
00351 << " by quadrature");
00352
00353 TEUCHOS_TEST_FOR_EXCEPTION(vectorCoeffs[resultIndices_[i]]->length()==0,
00354 std::logic_error,
00355 "zero-length coeff vector detected in "
00356 "quadrature integration branch of "
00357 "IntegralGroup::evaluate(). std::string value is ["
00358 << vectorCoeffs[resultIndices_[i]]->str()
00359 << "]");
00360
00361 Tabs tab3;
00362 SUNDANCE_MSG3(integrationVerb(),
00363 tab3 << "coefficients are " << vectorCoeffs[resultIndices_[i]]->str());
00364
00365 const double* const f = vectorCoeffs[resultIndices_[i]]->start();
00366 maxQuad->transform(JTrans, JVol, isLocalFlag, facetIndex, cellLIDs , f, A);
00367 }
00368 else if (curveQuad != 0)
00369 {
00370 SUNDANCE_MSG2(integrationVerb(),
00371 tab2 << "Integrating term group " << i
00372 << " by curve integral (quadrature by default) , result index: " << resultIndices_[i]);
00373
00374 double f_const = 0.0;
00375 if (constantCoeffs.size() > resultIndices_[i]){
00376 f_const = constantCoeffs[resultIndices_[i]];
00377 }
00378
00379 SUNDANCE_MSG2(integrationVerb(),
00380 tab2 << "Coefficient is " << f_const);
00381
00382
00383 if (vectorCoeffs.size() > resultIndices_[i]){
00384 Tabs tab3;
00385 double* const f = vectorCoeffs[resultIndices_[i]]->start();
00386 SUNDANCE_MSG3(integrationVerb(),
00387 tab3 << "coefficients are " << vectorCoeffs[resultIndices_[i]]->str());
00388 curveQuad->transform(JTrans, JVol, isLocalFlag, facetIndex, cellLIDs , f_const , f , A);
00389 } else{
00390 const double* f_null = 0;
00391 curveQuad->transform(JTrans, JVol, isLocalFlag, facetIndex, cellLIDs , f_const , f_null , A);
00392 }
00393
00394 }
00395 else
00396 {
00397 TEUCHOS_TEST_FOR_EXCEPT(1);
00398 }
00399
00400 SUNDANCE_MSG4(integrationVerb(),
00401 tab1 << "i=" << i << " integral values=");
00402 if (integrationVerb() >=4) writeTable(Out::os(), tab1, *A, 6);
00403 }
00404 SUNDANCE_MSG1(integrationVerb(), tab0 << "done integral group");
00405
00406 return true;
00407 }
00408
00409
00410 int IntegralGroup::findIntegrationVerb(const Array<RCP<ElementIntegral> >& integrals) const
00411 {
00412 int rtn = 0;
00413 for (int i=0; i<integrals.size(); i++)
00414 {
00415 rtn = std::max(rtn, integrals[i]->integrationVerb());
00416 }
00417 return rtn;
00418 }
00419
00420
00421 int IntegralGroup::findTransformVerb(const Array<RCP<ElementIntegral> >& integrals) const
00422 {
00423 int rtn = 0;
00424 for (int i=0; i<integrals.size(); i++)
00425 {
00426 rtn = std::max(rtn, integrals[i]->transformVerb());
00427 }
00428 return rtn;
00429 }
