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 "SundanceTrivialGrouper.hpp"
00032 #include "SundanceRefIntegral.hpp"
00033 #include "SundanceQuadratureIntegral.hpp"
00034 #include "SundanceReducedIntegral.hpp"
00035 #include "SundanceMaximalQuadratureIntegral.hpp"
00036 #include "SundanceCurveQuadratureIntegral.hpp"
00037 #include "SundanceEquationSet.hpp"
00038 #include "SundanceIntegralGroup.hpp"
00039 #include "SundanceBasisFamily.hpp"
00040 #include "SundanceSparsitySuperset.hpp"
00041 #include "SundanceQuadratureFamily.hpp"
00042 #include "SundanceReducedQuadrature.hpp"
00043 #include "SundanceMap.hpp"
00044 #include "SundanceOut.hpp"
00045 #include "PlayaTabs.hpp"
00046
00047 using namespace Sundance;
00048 using namespace Teuchos;
00049
00050
00051
00052 void TrivialGrouper::findGroups(const EquationSet& eqn,
00053 const CellType& maxCellType,
00054 int spatialDim,
00055 const CellType& cellType,
00056 int cellDim,
00057 const QuadratureFamily& quad,
00058 const RCP<SparsitySuperset>& sparsity,
00059 bool isInternalBdry,
00060 Array<RCP<IntegralGroup> >& groups,
00061 const ParametrizedCurve& globalCurve,
00062 const Mesh& mesh) const
00063 {
00064 Tabs tab(0);
00065
00066 SUNDANCE_MSG1(setupVerb(),
00067 tab << "in TrivialGrouper::findGroups(), num derivs = "
00068 << sparsity->numDerivs());
00069 SUNDANCE_MSG1(setupVerb(),
00070 tab << "cell type = " << cellType);
00071 SUNDANCE_MSG1(setupVerb(),
00072 tab << "sparsity = " << std::endl << *sparsity << std::endl);
00073
00074 const ReducedQuadrature* rq = dynamic_cast<const ReducedQuadrature*>(quad.ptr(
00075 ).get());
00076 bool useReducedQuad = (rq != 0);
00077 SUNDANCE_MSG1(setupVerb(), tab << "using reduced quadrature: "
00078 << useReducedQuad);
00079
00080
00081 int vecCount=0;
00082 int constCount=0;
00083
00084 bool isMaximal = cellType == maxCellType;
00085 bool useMaxIntegral = isMaximal;
00086
00087
00088
00089
00090 bool doGroups = true;
00091 if (!isMaximal) doGroups = false;
00092
00093
00094 typedef Sundance::Map<OrderedQuartet<int, BasisFamily, int, BasisFamily>, Array<RCP<ElementIntegral> > > twoFormMap;
00095 typedef Sundance::Map<OrderedTriple<int,int,BasisFamily>, Array<RCP<ElementIntegral> > > oneFormMap;
00096 Sundance::Map<OrderedQuartet<int, BasisFamily, int, BasisFamily>, Array<RCP<ElementIntegral> > > twoForms;
00097 Sundance::Map<OrderedQuartet<int, BasisFamily, int, BasisFamily>, Array<int> > twoFormResultIndices;
00098 Sundance::Map<OrderedTriple<int,int,BasisFamily>, Array<RCP<ElementIntegral> > > oneForms;
00099 Sundance::Map<OrderedTriple<int,int,BasisFamily>, Array<int> > oneFormResultIndices;
00100
00101 for (int i=0; i<sparsity->numDerivs(); i++)
00102 {
00103 const MultipleDeriv& d = sparsity->deriv(i);
00104 SUNDANCE_MSG3(setupVerb(),
00105 tab << "--------------------------------------------------");
00106 SUNDANCE_MSG3(setupVerb(),
00107 tab << "defining integration policy for " << d);
00108 SUNDANCE_MSG3(setupVerb(),
00109 tab << "--------------------------------------------------");
00110
00111 if (d.order()==0)
00112 {
00113 RCP<ElementIntegral> integral ;
00114 int resultIndex;
00115 if (sparsity->isConstant(i))
00116 {
00117 if (globalCurve.isCurveValid() && globalCurve.isCurveIntegral() && isMaximal )
00118 {
00119 integral = rcp(new CurveQuadratureIntegral( maxCellType, true ,
00120 quad, globalCurve , mesh , setupVerb() ) );
00121 }
00122 else
00123 {
00124 integral = rcp(new RefIntegral(spatialDim, maxCellType,
00125 cellDim, cellType, quad , isInternalBdry, globalCurve , mesh , setupVerb()));
00126 }
00127 resultIndex = constCount++;
00128 }
00129 else
00130 {
00131 if (useReducedQuad)
00132 {
00133 integral = rcp(new ReducedIntegral(spatialDim, maxCellType,
00134 cellDim, cellType, quad , isInternalBdry, globalCurve , mesh , setupVerb()));
00135 }
00136 else if (useMaxIntegral)
00137 {
00138 if (globalCurve.isCurveValid() && globalCurve.isCurveIntegral() && isMaximal )
00139 {
00140 integral = rcp(new CurveQuadratureIntegral(maxCellType, false ,
00141 quad, globalCurve , mesh , setupVerb()));
00142 }
00143 else
00144 {
00145 integral = rcp(new MaximalQuadratureIntegral(maxCellType,
00146 quad, globalCurve , mesh , setupVerb()));
00147 }
00148 }
00149 else
00150 {
00151 integral = rcp(new QuadratureIntegral(spatialDim, maxCellType,
00152 cellDim, cellType, quad, isInternalBdry, globalCurve , mesh ,
00153 setupVerb()));
00154 }
00155 resultIndex = vecCount++;
00156 }
00157 integral->setVerb(integrationVerb(), transformVerb());
00158 SUNDANCE_MSG3(setupVerb(), tab << "is zero-form");
00159 groups.append(rcp(new IntegralGroup(tuple(integral),
00160 tuple(resultIndex), setupVerb())));
00161 }
00162 else
00163 {
00164 BasisFamily testBasis;
00165 BasisFamily unkBasis;
00166 MultiIndex miTest;
00167 MultiIndex miUnk;
00168 int rawTestID = -1;
00169 int rawUnkID = -1;
00170 int rawParamID = -1;
00171 int testID = -1;
00172 int unkID = -1;
00173 int paramID = -1;
00174 int testBlock = -1;
00175 int unkBlock = -1;
00176 bool isOneForm;
00177 bool hasParam;
00178 extractWeakForm(eqn, d, testBasis, unkBasis, miTest, miUnk,
00179 rawTestID, rawUnkID,
00180 testID, unkID,
00181 testBlock, unkBlock,
00182 rawParamID, paramID,
00183 isOneForm, hasParam);
00184
00185 TEUCHOS_TEST_FOR_EXCEPT(hasParam && !isOneForm);
00186
00187
00188
00189
00190 int mvIndex = 0;
00191 if (hasParam) mvIndex = paramID;
00192
00193
00194
00195
00196
00197
00198
00199 bool transposeNeeded = false;
00200 if (!isOneForm && rawTestID!=rawUnkID
00201 && eqn.hasVarID(rawUnkID) && eqn.hasUnkID(rawTestID))
00202 {
00203 transposeNeeded = true;
00204 }
00205
00206
00207 if (isOneForm)
00208 {
00209 SUNDANCE_MSG3(setupVerb(), tab << "is one-form");
00210 }
00211 else
00212 {
00213 SUNDANCE_MSG3(setupVerb(), tab << "is two-form");
00214 }
00215
00216
00217 if (hasParam)
00218 {
00219 SUNDANCE_MSG3(setupVerb(), tab << "is a parametric derivative");
00220 }
00221 else
00222 {
00223 SUNDANCE_MSG3(setupVerb(), tab << "is not a parametric derivative");
00224 }
00225
00226 SUNDANCE_MSG3(setupVerb(),
00227 tab << "test ID: " << testID << " block=" << testBlock);
00228
00229 if (!isOneForm)
00230 {
00231 SUNDANCE_MSG3(setupVerb(), tab << "unk funcID: " << unkID << " block=" << unkBlock);
00232 }
00233
00234 if (hasParam)
00235 {
00236 SUNDANCE_MSG3(setupVerb(), tab << "paramID=" << paramID);
00237 }
00238
00239 SUNDANCE_MSG3(setupVerb(), tab << "deriv = " << d);
00240 if (sparsity->isConstant(i))
00241 {
00242 SUNDANCE_MSG3(setupVerb(), tab << "coeff is constant");
00243 }
00244 else
00245 {
00246 SUNDANCE_MSG3(setupVerb(), tab << "coeff is non-constant");
00247 }
00248
00249 RCP<ElementIntegral> integral;
00250 RCP<ElementIntegral> transposedIntegral;
00251 int resultIndex;
00252 if (sparsity->isConstant(i))
00253 {
00254 if (isOneForm)
00255 {
00256 int alpha=0;
00257 if (miTest.order()==1)
00258 {
00259 alpha = miTest.firstOrderDirection();
00260 }
00261
00262 if (globalCurve.isCurveValid() && globalCurve.isCurveIntegral() && isMaximal )
00263 {
00264 integral = rcp(new CurveQuadratureIntegral(maxCellType, true ,
00265 testBasis, alpha,
00266 miTest.order(), quad, globalCurve , mesh ,setupVerb()));
00267 }
00268 else
00269 {
00270 integral = rcp(new RefIntegral(spatialDim, maxCellType,
00271 cellDim, cellType,
00272 testBasis, alpha,
00273 miTest.order(), quad , isInternalBdry, globalCurve , mesh ,setupVerb()));
00274 }
00275 }
00276 else
00277 {
00278 int alpha=0;
00279 int beta=0;
00280 if (miTest.order()==1)
00281 {
00282 alpha = miTest.firstOrderDirection();
00283 }
00284 if (miUnk.order()==1)
00285 {
00286 beta = miUnk.firstOrderDirection();
00287 }
00288
00289 if (globalCurve.isCurveValid() && globalCurve.isCurveIntegral() && isMaximal )
00290 {
00291 integral = rcp(new CurveQuadratureIntegral(maxCellType, true ,
00292 testBasis, alpha,
00293 miTest.order(),
00294 unkBasis, beta,
00295 miUnk.order(), quad, globalCurve , mesh ,setupVerb()));
00296 }
00297 else
00298 {
00299 integral = rcp(new RefIntegral(spatialDim, maxCellType,
00300 cellDim, cellType,
00301 testBasis, alpha, miTest.order(),
00302 unkBasis, beta, miUnk.order(), quad , isInternalBdry, globalCurve , mesh ,setupVerb()));
00303 }
00304 if (transposeNeeded)
00305 {
00306 if (globalCurve.isCurveValid() && globalCurve.isCurveIntegral() && isMaximal )
00307 {
00308 transposedIntegral = rcp(new CurveQuadratureIntegral(maxCellType, true ,
00309 unkBasis, beta,
00310 miUnk.order(),
00311 testBasis, alpha,
00312 miTest.order(),
00313 quad, globalCurve , mesh ,setupVerb()));
00314 }
00315 else
00316 {
00317 transposedIntegral = rcp(new RefIntegral(spatialDim, maxCellType,
00318 cellDim, cellType,
00319 unkBasis, beta,
00320 miUnk.order(),
00321 testBasis, alpha,
00322 miTest.order(), quad , isInternalBdry, globalCurve , mesh ,setupVerb()));
00323 }
00324 }
00325 }
00326 resultIndex = constCount++;
00327 }
00328 else
00329 {
00330 if (isOneForm)
00331 {
00332 int alpha=0;
00333 if (miTest.order()==1)
00334 {
00335 alpha = miTest.firstOrderDirection();
00336 }
00337 if (useReducedQuad)
00338 {
00339 integral = rcp(new ReducedIntegral(spatialDim, maxCellType,
00340 cellDim, cellType,
00341 testBasis, alpha,
00342 miTest.order(), quad , isInternalBdry, globalCurve , mesh ,setupVerb()));
00343 }
00344 else if (useMaxIntegral)
00345 {
00346 if ( globalCurve.isCurveValid() && globalCurve.isCurveIntegral() )
00347 {
00348 integral = rcp(new CurveQuadratureIntegral(maxCellType, false ,
00349 testBasis, alpha,
00350 miTest.order(), quad, globalCurve , mesh ,setupVerb()));
00351 }
00352 else
00353 {
00354 integral = rcp(new MaximalQuadratureIntegral(maxCellType,
00355 testBasis, alpha,
00356 miTest.order(), quad, globalCurve , mesh ,setupVerb()));
00357 }
00358 }
00359 else
00360 {
00361 integral = rcp(new QuadratureIntegral(spatialDim, maxCellType,
00362 cellDim, cellType,
00363 testBasis, alpha,
00364 miTest.order(), quad, isInternalBdry, globalCurve , mesh ,setupVerb()));
00365 }
00366 }
00367 else
00368 {
00369 int alpha=0;
00370 int beta=0;
00371 if (miTest.order()==1)
00372 {
00373 alpha = miTest.firstOrderDirection();
00374 }
00375 if (miUnk.order()==1)
00376 {
00377 beta = miUnk.firstOrderDirection();
00378 }
00379 if (useReducedQuad)
00380 {
00381 integral = rcp(new ReducedIntegral(spatialDim, maxCellType,
00382 cellDim, cellType,
00383 testBasis, alpha, miTest.order(),
00384 unkBasis, beta, miUnk.order(), quad , isInternalBdry, globalCurve , mesh ,setupVerb()));
00385 if (transposeNeeded)
00386 {
00387 transposedIntegral = rcp(new ReducedIntegral(spatialDim, maxCellType,
00388 cellDim, cellType,
00389 unkBasis, beta,
00390 miUnk.order(),
00391 testBasis, alpha,
00392 miTest.order(), quad , isInternalBdry, globalCurve , mesh ,setupVerb()));
00393 }
00394 }
00395 else if (useMaxIntegral)
00396 {
00397 if ( globalCurve.isCurveValid() && globalCurve.isCurveIntegral() )
00398 {
00399 integral = rcp(new CurveQuadratureIntegral(maxCellType, false ,
00400 testBasis, alpha,
00401 miTest.order(),
00402 unkBasis, beta,
00403 miUnk.order(), quad, globalCurve , mesh ,setupVerb()));
00404 }
00405 else
00406 {
00407 integral = rcp(new MaximalQuadratureIntegral(maxCellType,
00408 testBasis, alpha,
00409 miTest.order(),
00410 unkBasis, beta,
00411 miUnk.order(), quad, globalCurve , mesh ,setupVerb()));
00412 }
00413 if (transposeNeeded)
00414 {
00415 if ( globalCurve.isCurveValid() && globalCurve.isCurveIntegral() )
00416 {
00417 transposedIntegral = rcp(new CurveQuadratureIntegral(maxCellType, false ,
00418 unkBasis, beta, miUnk.order(),
00419 testBasis, alpha, miTest.order(), quad,
00420 globalCurve , mesh ,setupVerb()));
00421 }
00422 else
00423 {
00424 transposedIntegral = rcp(new MaximalQuadratureIntegral(maxCellType,
00425 unkBasis, beta, miUnk.order(),
00426 testBasis, alpha, miTest.order(), quad,
00427 globalCurve , mesh ,setupVerb()));
00428 }
00429 }
00430 }
00431 else
00432 {
00433 integral = rcp(new QuadratureIntegral(spatialDim, maxCellType,
00434 cellDim, cellType,
00435 testBasis, alpha,
00436 miTest.order(),
00437 unkBasis, beta,
00438 miUnk.order(), quad, isInternalBdry, globalCurve , mesh ,setupVerb()));
00439 if (transposeNeeded)
00440 {
00441 transposedIntegral = rcp(new QuadratureIntegral(spatialDim, maxCellType,
00442 cellDim, cellType,
00443 unkBasis, beta, miUnk.order(),
00444 testBasis, alpha, miTest.order(), quad,
00445 isInternalBdry, globalCurve , mesh ,setupVerb()));
00446 }
00447 }
00448 }
00449 resultIndex = vecCount++;
00450 }
00451
00452
00453 integral->setVerb(integrationVerb(), transformVerb());
00454 if (transposeNeeded)
00455 {
00456 transposedIntegral->setVerb(integrationVerb(), transformVerb());
00457 }
00458
00459
00460 if (isOneForm)
00461 {
00462 if (doGroups)
00463 {
00464 OrderedTriple<int,int,BasisFamily> testKey(rawTestID, mvIndex, testBasis);
00465 if (!oneForms.containsKey(testKey))
00466 {
00467 oneForms.put(testKey, tuple(integral));
00468 oneFormResultIndices.put(testKey, tuple(resultIndex));
00469 }
00470 else
00471 {
00472 oneForms[testKey].append(integral);
00473 oneFormResultIndices[testKey].append(resultIndex);
00474 }
00475 }
00476 else
00477 {
00478 groups.append(rcp(new IntegralGroup(tuple(testID), tuple(testBlock),
00479 tuple(mvIndex),
00480 tuple(integral),
00481 tuple(resultIndex), tuple(d), setupVerb())));
00482 }
00483 }
00484 else
00485 {
00486 if (!doGroups)
00487 {
00488 groups.append(rcp(new IntegralGroup(tuple(testID), tuple(testBlock),
00489 tuple(unkID), tuple(unkBlock),
00490 tuple(integral),
00491 tuple(resultIndex), tuple(d), setupVerb())));
00492 if (transposeNeeded)
00493 {
00494 groups.append(rcp(new IntegralGroup(tuple(unkID), tuple(unkBlock),
00495 tuple(testID), tuple(testBlock),
00496 tuple(transposedIntegral),
00497 tuple(resultIndex), tuple(d), setupVerb())));
00498 }
00499 }
00500 else
00501 {
00502 Tabs tab3;
00503 OrderedQuartet<int, BasisFamily, int, BasisFamily> testUnkKey(rawTestID, testBasis, rawUnkID, unkBasis);
00504
00505
00506 SUNDANCE_MSG2(setupVerb(), tab3 << "key=" << testUnkKey);
00507 if (!twoForms.containsKey(testUnkKey))
00508 {
00509 Tabs tab4;
00510 SUNDANCE_MSG2(setupVerb(), tab4 << "key not found");
00511 twoForms.put(testUnkKey, tuple(integral));
00512 twoFormResultIndices.put(testUnkKey, tuple(resultIndex));
00513 }
00514 else
00515 {
00516 Tabs tab4;
00517 SUNDANCE_MSG2(setupVerb(), tab4 << "key found");
00518 twoForms[testUnkKey].append(integral);
00519 twoFormResultIndices[testUnkKey].append(resultIndex);
00520 }
00521 if (transposeNeeded)
00522 {
00523 OrderedQuartet<int, BasisFamily, int, BasisFamily> unkTestKey(rawUnkID, unkBasis, rawTestID, testBasis);
00524
00525 if (!twoForms.containsKey(unkTestKey))
00526 {
00527 Tabs tab4;
00528 SUNDANCE_MSG2(setupVerb(), tab4 << "key not found");
00529 twoForms.put(unkTestKey, tuple(transposedIntegral));
00530 twoFormResultIndices.put(unkTestKey, tuple(resultIndex));
00531 }
00532 else
00533 {
00534 Tabs tab4;
00535 SUNDANCE_MSG2(setupVerb(), tab4 << "key found");
00536 twoForms[unkTestKey].append(transposedIntegral);
00537 twoFormResultIndices[unkTestKey].append(resultIndex);
00538 }
00539 }
00540 }
00541 }
00542 }
00543 }
00544
00545 if (doGroups)
00546 {
00547 Tabs tab;
00548 SUNDANCE_MSG2(setupVerb(), tab << "creating integral groups");
00549 for (twoFormMap::const_iterator i=twoForms.begin(); i!=twoForms.end(); i++)
00550 {
00551 Tabs tab3;
00552 SUNDANCE_MSG2(setupVerb(), tab3 << "integral group number="
00553 << groups.size());
00554 int rawTestID = i->first.a();
00555 BasisFamily testBasis = i->first.b();
00556 int rawUnkID = i->first.c();
00557 BasisFamily unkBasis = i->first.d();
00558 int testID = eqn.reducedVarID(rawTestID);
00559 int unkID = eqn.reducedUnkID(rawUnkID);
00560 int testBlock = eqn.blockForVarID(rawTestID);
00561 int unkBlock = eqn.blockForUnkID(rawUnkID);
00562 const Array<RCP<ElementIntegral> >& integrals = i->second;
00563 const Array<int>& resultIndices
00564 = twoFormResultIndices.get(i->first);
00565 SUNDANCE_MSG2(setupVerb(), tab3 << "creating two-form integral group" << std::endl
00566 << tab3 << "testID=" << rawTestID << std::endl
00567 << tab3 << "unkID=" << rawUnkID << std::endl
00568 << tab3 << "testBlock=" << testBlock << std::endl
00569 << tab3 << "unkBlock=" << unkBlock << std::endl
00570 << tab3 << "testBasis=" << testBasis << std::endl
00571 << tab3 << "unkBasis=" << unkBasis << std::endl
00572 << tab3 << "resultIndices=" << resultIndices);
00573 Array<MultipleDeriv> grpDerivs;
00574 for (int j=0; j<resultIndices.size(); j++)
00575 {
00576 MultipleDeriv d = sparsity->deriv(resultIndices[j]);
00577 SUNDANCE_MSG2(setupVerb(), tab3 << "deriv " << j << " "
00578 << d);
00579 grpDerivs.append(d);
00580 }
00581 groups.append(rcp(new IntegralGroup(tuple(testID), tuple(testBlock),
00582 tuple(unkID), tuple(unkBlock),
00583 integrals, resultIndices, grpDerivs, setupVerb())));
00584 }
00585
00586 for (oneFormMap::const_iterator i=oneForms.begin(); i!=oneForms.end(); i++)
00587 {
00588 Tabs tab3;
00589 SUNDANCE_MSG2(setupVerb(), tab3 << "integral group number="
00590 << groups.size());
00591 int rawTestID = i->first.a();
00592 int mvIndex = i->first.b();
00593 int testID = eqn.reducedVarID(rawTestID);
00594 int testBlock = eqn.blockForVarID(rawTestID);
00595 const Array<RCP<ElementIntegral> >& integrals = i->second;
00596 const Array<int>& resultIndices
00597 = oneFormResultIndices.get(i->first);
00598 SUNDANCE_MSG2(setupVerb(), tab3 << "creating one-form integral group" << std::endl
00599 << tab3 << "testID=" << testID << std::endl
00600 << tab3 << "resultIndices=" << resultIndices);
00601 Array<MultipleDeriv> grpDerivs;
00602 for (int j=0; j<resultIndices.size(); j++)
00603 {
00604 MultipleDeriv d = sparsity->deriv(resultIndices[j]);
00605 SUNDANCE_MSG2(setupVerb(), tab3 << "deriv " << j << " "
00606 << d);
00607 grpDerivs.append(d);
00608 }
00609 groups.append(rcp(new IntegralGroup(tuple(testID), tuple(testBlock),
00610 tuple(mvIndex),
00611 integrals, resultIndices, grpDerivs, setupVerb())));
00612 }
00613 }
00614
00615
00616 }
00617
