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