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