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 "SundanceFeketeBrickQuadrature.hpp"
00032 #include "SundanceFeketeQuadQuadrature.hpp"
00033 #include "SundanceFeketeQuadrature.hpp"
00034 #include "SundanceFeketeTriangleQuadrature.hpp"
00035 #include "SundanceGaussLobatto1D.hpp"
00036 #include "SundanceOut.hpp"
00037 #include "PlayaTabs.hpp"
00038 #include <stack>
00039
00040 using namespace Sundance;
00041 using namespace Teuchos;
00042
00043
00044 extern "C"
00045 {
00046
00047 void dgemv_(char *trans, int *m, int *n, double *alpha, double *a, int *lda,
00048 double *x, int *incx, double *beta, double *y, int *incy);
00049 }
00050
00051
00052 FeketeQuadrature::FeketeQuadrature(int order) :
00053 QuadratureFamilyBase(order)
00054 {
00055 _hasBasisCoeffs = false;
00056 }
00057
00058 XMLObject FeketeQuadrature::toXML() const
00059 {
00060 XMLObject rtn("FeketeQuadrature");
00061 rtn.addAttribute("order", Teuchos::toString(order()));
00062 return rtn;
00063 }
00064
00065 void FeketeQuadrature::getLineRule(Array<Point>& quadPoints,
00066 Array<double>& quadWeights) const
00067 {
00068 int p = order() + 3;
00069 p = p + (p % 2);
00070 int n = p / 2;
00071
00072 quadPoints.resize(n);
00073 quadWeights.resize(n);
00074
00075 GaussLobatto1D q1(n, 0.0, 1.0);
00076
00077 for (int i = 0; i < n; i++)
00078 {
00079 quadWeights[i] = q1.weights()[i];
00080 quadPoints[i] = Point(q1.nodes()[i]);
00081 }
00082 }
00083
00084 void FeketeQuadrature::getTriangleRule(Array<Point>& quadPoints,
00085 Array<double>& quadWeights) const
00086 {
00087 Array<double> x;
00088 Array<double> y;
00089 Array<double> w;
00090
00091 FeketeTriangleQuadrature::getPoints(order(), w, x, y);
00092 quadPoints.resize(w.length());
00093 quadWeights.resize(w.length());
00094 for (int i = 0; i < w.length(); i++)
00095 {
00096 quadWeights[i] = 0.5 * w[i];
00097 quadPoints[i] = Point(x[i], y[i]);
00098 }
00099 }
00100
00101 void FeketeQuadrature::getQuadRule(Array<Point>& quadPoints,
00102 Array<double>& quadWeights) const
00103 {
00104 Array<double> x;
00105 Array<double> y;
00106 Array<double> w;
00107
00108 FeketeQuadQuadrature::getPoints(order(), w, x, y);
00109 quadPoints.resize(w.length());
00110 quadWeights.resize(w.length());
00111 for (int i = 0; i < w.length(); i++)
00112 {
00113 quadWeights[i] = w[i];
00114 quadPoints[i] = Point(x[i], y[i]);
00115 }
00116 }
00117
00118 void FeketeQuadrature::getBrickRule(Array<Point>& quadPoints,
00119 Array<double>& quadWeights) const
00120 {
00121 Array<double> x;
00122 Array<double> y;
00123 Array<double> z;
00124 Array<double> w;
00125
00126 FeketeBrickQuadrature::getPoints(order(), w, x, y, z);
00127 quadPoints.resize(w.length());
00128 quadWeights.resize(w.length());
00129 for (int i = 0; i < w.length(); i++)
00130 {
00131 quadWeights[i] = w[i];
00132 quadPoints[i] = Point(x[i], y[i], z[i]);
00133 }
00134 }
00135
00136 void FeketeQuadrature::getAdaptedWeights(const CellType& cellType, int cellDim,
00137 int cellLID, int facetIndex, const Mesh& mesh, const ParametrizedCurve&
00138 globalCurve, Array<Point>& quadPoints, Array<double>& quadWeights,
00139 bool& weightsChanged) const
00140 {
00141
00142 Tabs tabs;
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 CellType maxCellType = mesh.cellType(mesh.spatialDim());
00155
00156 switch (maxCellType)
00157 {
00158
00159 case TriangleCell:
00160 {
00161 switch (cellType)
00162 {
00163 case TriangleCell:
00164 {
00165 getAdaptedTriangleWeights(cellLID, mesh, globalCurve, quadPoints,
00166 quadWeights, weightsChanged);
00167 break;
00168 }
00169 default:
00170 #ifndef TRILINOS_7
00171 SUNDANCE_ERROR ("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for submaximal cell " << maxCellType << " in a triangle mesh")
00172 ;
00173 #else
00174 SUNDANCE_ERROR7("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for submaximal cell " << maxCellType << " in a triangle mesh");
00175 #endif
00176 }
00177 break;
00178 }
00179 case QuadCell:
00180 {
00181 switch(cellType)
00182 {
00183 case QuadCell:
00184 {
00185 getAdaptedQuadWeights(cellLID, mesh, globalCurve, quadPoints,
00186 quadWeights, weightsChanged);
00187 break;
00188 }
00189 default:
00190 #ifndef TRILINOS_7
00191 SUNDANCE_ERROR ("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for submaximal cell " << maxCellType << " in a triangle mesh")
00192 ;
00193 #else
00194 SUNDANCE_ERROR7("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for submaximal cell " << maxCellType << " in a triangle mesh");
00195 #endif
00196 }
00197 break;
00198 }
00199 default:
00200 #ifndef TRILINOS_7
00201 SUNDANCE_ERROR("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for cell type " << cellType)
00202 ;
00203 #else
00204 SUNDANCE_ERROR7("getAdaptedWeights rule for ACI(Adaptive Cell Integration) not available for cell type " << cellType);
00205 #endif
00206 }
00207 }
00208
00209 void FeketeQuadrature::getAdaptedTriangleWeights(int cellLID, const Mesh& mesh,
00210 const ParametrizedCurve& globalCurve, Array<Point>& quadPoints, Array<
00211 double>& quadWeights, bool& weightsChanged) const
00212 {
00213 int verb = 0;
00214 Tabs tabs;
00215
00216
00217 double tol = 1.e-8;
00218 weightsChanged = true;
00219
00220
00221 Array<int> cellLIDs(1);
00222 cellLIDs[0] = cellLID;
00223
00224
00225 Array<double> ip01, ip02, ip12;
00226 int nr_ip01, nr_ip02, nr_ip12;
00227
00228 Array<double> originalWeights = quadWeights;
00229 ParametrizedCurve noCurve = new DummyParametrizedCurve();
00230 Array<double> integrals(originalWeights.size());
00231 for (int i = 0; i < integrals.size(); i++)
00232 integrals[i] = 0.0;
00233
00234
00235 std::stack<Point> s;
00236 std::stack<double> sa;
00237 s.push(Point(0.0, 0.0));
00238 s.push(Point(1.0, 0.0));
00239 s.push(Point(0.0, 1.0));
00240 sa.push(0.5);
00241 Array<Point> nodes;
00242 Array<Point> nodesref(3);
00243
00244
00245 unsigned int nTris = 0;
00246
00247 while (s.size() >= 3)
00248 {
00249
00250 ++nTris;
00251
00252 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedTriangleWeights Current number of triangles on stack: " << sa.size());
00253
00254
00255 bool refine = false;
00256 Array<double> results1(originalWeights.size());
00257 Array<double> results2(originalWeights.size());
00258
00259
00260 double area = sa.top();
00261 sa.pop();
00262 for (int i = 2; i >= 0; i--)
00263 {
00264 nodesref[i] = s.top();
00265 s.pop();
00266 }
00267 mesh.pushForward(2, cellLIDs, nodesref, nodes);
00268
00269
00270 globalCurve.returnIntersect(nodes[0], nodes[1], nr_ip01, ip01);
00271 globalCurve.returnIntersect(nodes[0], nodes[2], nr_ip02, ip02);
00272 globalCurve.returnIntersect(nodes[1], nodes[2], nr_ip12, ip12);
00273
00274
00275 Point x, xref;
00276 Point vec1, vec1ref;
00277 Point vec2, vec2ref;
00278 bool xInOmega = false;
00279
00280
00281 if ((nr_ip01 + nr_ip02 + nr_ip12 == 0) || (nr_ip01 == 0 && nr_ip02 == 0
00282 && nr_ip12 == 1) || (nr_ip01 == 1 && nr_ip02 == 1 && nr_ip12
00283 == 0))
00284 {
00285 x = nodes[0];
00286 vec1 = nodes[2] - nodes[0];
00287 vec2 = nodes[1] - nodes[0];
00288 xref = nodesref[0];
00289 vec1ref = nodesref[2] - nodesref[0];
00290 vec2ref = nodesref[1] - nodesref[0];
00291 }
00292 else if ((nr_ip01 == 0 && nr_ip02 == 1 && nr_ip12 == 0) || (nr_ip12
00293 == 1 && nr_ip01 == 1 && nr_ip02 == 0))
00294 {
00295 x = nodes[1];
00296 vec1 = nodes[0] - nodes[1];
00297 vec2 = nodes[2] - nodes[1];
00298 xref = nodesref[1];
00299 vec1ref = nodesref[0] - nodesref[1];
00300 vec2ref = nodesref[2] - nodesref[1];
00301 }
00302 else if ((nr_ip01 == 1 && nr_ip02 == 0 && nr_ip12 == 0) || (nr_ip02
00303 == 1 && nr_ip12 == 1 && nr_ip01 == 0))
00304 {
00305 x = nodes[2];
00306 vec1 = nodes[1] - nodes[2];
00307 vec2 = nodes[0] - nodes[2];
00308 xref = nodesref[2];
00309 vec1ref = nodesref[1] - nodesref[2];
00310 vec2ref = nodesref[0] - nodesref[2];
00311 }
00312 else
00313 {
00314 refine = true;
00315 }
00316
00317
00318 if (!refine)
00319 {
00320 if (globalCurve.curveEquation(x) < 0)
00321 xInOmega = true;
00322 integrateRegion(TriangleCell, 2, order() + 1, x, xref, vec1, vec2,
00323 vec1ref, vec2ref, globalCurve, results1);
00324 if (results1.size() == 0)
00325 {
00326 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedTriangleWeights Unknown error in integration. Refine.")
00327 refine = true;
00328 }
00329 else
00330 {
00331 integrateRegion(TriangleCell, 2, order() + 2, x, xref, vec1,
00332 vec2, vec1ref, vec2ref, globalCurve, results2);
00333
00334
00335 for (int i = 0; i < results1.size(); i++)
00336 {
00337 if (::fabs(results2[i] - results1[i]) > sqrt(area * tol))
00338 {
00339 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedTriangleWeights Insufficient accuracy. Refine.")
00340 refine = true;
00341 break;
00342 }
00343 }
00344 }
00345 }
00346
00347
00348 if (refine)
00349 {
00350
00351
00352 s.push(nodesref[0]);
00353 s.push((nodesref[0] + nodesref[1]) / 2);
00354 s.push((nodesref[0] + nodesref[2]) / 2);
00355 sa.push(area / 4);
00356 s.push((nodesref[0] + nodesref[1]) / 2);
00357 s.push(nodesref[1]);
00358 s.push((nodesref[1] + nodesref[2]) / 2);
00359 sa.push(area / 4);
00360 s.push((nodesref[0] + nodesref[2]) / 2);
00361 s.push((nodesref[1] + nodesref[2]) / 2);
00362 s.push(nodesref[2]);
00363 sa.push(area / 4);
00364 s.push((nodesref[0] + nodesref[2]) / 2);
00365 s.push((nodesref[0] + nodesref[1]) / 2);
00366 s.push((nodesref[1] + nodesref[2]) / 2);
00367 sa.push(area / 4);
00368 continue;
00369 }
00370
00371
00372 if (xInOmega)
00373 {
00374 for (int i = 0; i < integrals.size(); i++)
00375 integrals[i] += results2[i];
00376 }
00377 else
00378 {
00379 integrateRegion(TriangleCell, 2, order() + 1, x, xref, vec1, vec2,
00380 vec1ref, vec2ref, noCurve, results1);
00381 for (int i = 0; i < integrals.size(); i++)
00382 integrals[i] += (results1[i] - results2[i]);
00383 }
00384 }
00385
00386 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedTriangleWeights Calculations in cell " << cellLID << " completed. " << nTris << " triangle(s) used.");
00387
00388 Array<double> alphas;
00389 globalCurve.getIntegrationParams(alphas);
00390 for (int i = 0; i < quadWeights.size(); i++)
00391 quadWeights[i] = alphas[1] * integrals[i] + alphas[0]
00392 * (originalWeights[i] - integrals[i]);
00393 mesh.setSpecialWeight(2, cellLID, quadWeights);
00394 }
00395
00396 void FeketeQuadrature::getAdaptedQuadWeights(int cellLID, const Mesh& mesh,
00397 const ParametrizedCurve& globalCurve, Array<Point>& quadPoints, Array<
00398 double>& quadWeights, bool& weightsChanged) const
00399 {
00400
00401 int verb = 0;
00402 Tabs tabs;
00403
00404
00405 double tol = 1.e-8;
00406
00407
00408 weightsChanged = true;
00409 Array<double> originalWeights = quadWeights;
00410
00411 ParametrizedCurve noCurve = new DummyParametrizedCurve();
00412 Array<double> integrals(originalWeights.size());
00413 for (int i = 0; i < integrals.size(); i++)
00414 integrals[i] = 0.0;
00415
00416
00417 Array<int> cellLIDs(1);
00418 cellLIDs[0] = cellLID;
00419 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedQuadWeights: Current cell LID " << cellLID)
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430 Array<int> nIntEdge(4);
00431 Array<Array<double> > IntEdgePars(4);
00432 int edgeIndex[4][2] = { {0,1} , {0,2} , {1,3} , {2,3} };
00433
00434
00435 std::stack<Point> s;
00436 std::stack<double> sa;
00437 s.push(Point(0.0, 0.0));
00438 s.push(Point(1.0, 0.0));
00439 s.push(Point(0.0, 1.0));
00440 s.push(Point(1.0, 1.0));
00441 sa.push(1.);
00442 Array<Point> nodes;
00443 Array<Point> nodesref(4);
00444
00445
00446 unsigned int nQuads = 0;
00447
00448 while (s.size() >= 4)
00449 {
00450
00451 ++nQuads;
00452 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedQuadWeights: Current number of quads on stack: " << sa.size())
00453
00454
00455 bool refine = false;
00456 bool xInOmega = false;
00457 Array<double> results1(originalWeights.size());
00458 Array<double> results2(originalWeights.size());
00459
00460
00461 double area = sa.top();
00462 sa.pop();
00463 for (int i = 3; i >= 0; i--)
00464 {
00465 nodesref[i] = s.top();
00466 s.pop();
00467 }
00468 mesh.pushForward(2, cellLIDs, nodesref, nodes);
00469
00470
00471 for (int jj = 0 ; jj < 4 ; jj++ ) {
00472 globalCurve.returnIntersect(nodes[edgeIndex[jj][0]], nodes[edgeIndex[jj][1]], nIntEdge[jj], IntEdgePars[jj] );
00473
00474
00475
00476
00477
00478 }
00479
00480
00481 Point x, xref;
00482 Point vec1, vec1ref;
00483 Point vec2, vec2ref;
00484
00485
00486 if ( nIntEdge[0] == 0 )
00487 {
00488 x = nodes[0];
00489 vec1 = nodes[1] - nodes[0];
00490 vec2 = nodes[2] - nodes[0];
00491 xref = nodesref[0];
00492 vec1ref = nodesref[1] - nodesref[0];
00493 vec2ref = nodesref[2] - nodesref[0];
00494 }
00495 else if ( nIntEdge[1] == 0 )
00496 {
00497 x = nodes[2];
00498 vec1 = nodes[0] - nodes[2];
00499 vec2 = nodes[3] - nodes[2];
00500 xref = nodesref[2];
00501 vec1ref = nodesref[0] - nodesref[2];
00502 vec2ref = nodesref[3] - nodesref[2];
00503 }
00504 else if ( nIntEdge[2] == 0 )
00505 {
00506 x = nodes[1];
00507 vec1 = nodes[3] - nodes[1];
00508 vec2 = nodes[0] - nodes[1];
00509 xref = nodesref[1];
00510 vec1ref = nodesref[3] - nodesref[1];
00511 vec2ref = nodesref[0] - nodesref[1];
00512 }
00513 else
00514 {
00515 refine = true;
00516 }
00517
00518
00519 if (!refine)
00520 {
00521 if (globalCurve.curveEquation(x) < 0)
00522 xInOmega = true;
00523 integrateRegion(QuadCell, 2, order() + 1, x, xref, vec1, vec2,
00524 vec1ref, vec2ref, globalCurve, results1);
00525 if (results1.size() == 0)
00526 {
00527 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedQuadWeights Unknown error in integration. Refine.")
00528 refine = true;
00529 }
00530 else
00531 {
00532 integrateRegion(QuadCell, 2, order() + 2, x, xref, vec1,
00533 vec2, vec1ref, vec2ref, globalCurve, results2);
00534
00535
00536 double reltol = ::sqrt(area * tol);
00537 for (int i = 0; i < results1.size(); i++)
00538 {
00539 if (::fabs(results2[i] - results1[i]) > reltol)
00540 {
00541 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedQuadWeights Insufficient accuracy. Refine.")
00542 refine = true;
00543 break;
00544 }
00545 }
00546 }
00547 }
00548
00549
00550 if (refine)
00551 {
00552
00553
00554 s.push(nodesref[0]);
00555 s.push((nodesref[0] + nodesref[1]) / 2);
00556 s.push((nodesref[0] + nodesref[2]) / 2);
00557 s.push((nodesref[0] + nodesref[3]) / 2);
00558 sa.push(area / 4);
00559 s.push((nodesref[0] + nodesref[1]) / 2);
00560 s.push(nodesref[1]);
00561 s.push((nodesref[1] + nodesref[2]) / 2);
00562 s.push((nodesref[1] + nodesref[3]) / 2);
00563 sa.push(area / 4);
00564 s.push((nodesref[0] + nodesref[2]) / 2);
00565 s.push((nodesref[1] + nodesref[2]) / 2);
00566 s.push(nodesref[2]);
00567 s.push((nodesref[2] + nodesref[3]) / 2);
00568 sa.push(area / 4);
00569 s.push((nodesref[0] + nodesref[3]) / 2);
00570 s.push((nodesref[1] + nodesref[3]) / 2);
00571 s.push((nodesref[2] + nodesref[3]) / 2);
00572 s.push(nodesref[3]);
00573 sa.push(area / 4);
00574 continue;
00575 }
00576
00577
00578 if (xInOmega)
00579 {
00580 for (int i = 0; i < integrals.size(); i++)
00581 integrals[i] += results2[i];
00582 }
00583 else
00584 {
00585 integrateRegion(QuadCell, 2, order() + 1, x, xref, vec1, vec2,
00586 vec1ref, vec2ref, noCurve, results1);
00587 for (int i = 0; i < integrals.size(); i++)
00588 integrals[i] += (results1[i] - results2[i]);
00589 }
00590 }
00591
00592 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::getAdaptedQuadWeights Calculations in cell " << cellLID << " completed. " << nQuads << " quad(s) used.");
00593
00594 Array<double> alphas;
00595 globalCurve.getIntegrationParams(alphas);
00596 for (int i = 0; i < quadWeights.size(); i++)
00597 quadWeights[i] = alphas[1] * integrals[i] + alphas[0]
00598 * (originalWeights[i] - integrals[i]);
00599 if (verb > 1 ) {
00600 writeTable(Out::os(), tabs, quadWeights, quadWeights.size());
00601 }
00602 mesh.setSpecialWeight(2, cellLID, quadWeights);
00603 }
00604
00605 void FeketeQuadrature::integrateRegion(const CellType& cellType,
00606 const int cellDim, const int innerOrder, const Point& x,
00607 const Point& xref, const Point& vec1, const Point& vec2,
00608 const Point& vec1ref, const Point& vec2ref,
00609 const ParametrizedCurve& curve, Array<double>& integrals) const
00610 {
00611 int verb = 0;
00612 Tabs tabs;
00613
00614 switch (cellType)
00615 {
00616 case TriangleCell:
00617 {
00618
00619 integrals.resize((order() + 1) * (order() + 2) / 2);
00620 for (int i = 0; i < integrals.size(); i++)
00621 integrals[i] = 0.0;
00622
00623
00624 int nrInnerQuadPts = innerOrder + 3;
00625 nrInnerQuadPts = (nrInnerQuadPts + nrInnerQuadPts % 2) / 2;
00626 GaussLobatto1D innerQuad(nrInnerQuadPts, 0.0, 1.0);
00627
00628 Array<double> innerQuadPts(nrInnerQuadPts);
00629 Array<double> innerQuadWgts(nrInnerQuadPts);
00630 innerQuadPts = innerQuad.nodes();
00631 innerQuadWgts = innerQuad.weights();
00632
00633
00634
00635 for (int edgePos = 0; edgePos < nrInnerQuadPts; edgePos++)
00636 {
00637
00638 Point ray = innerQuadPts[edgePos] * vec1 + (1.0
00639 - innerQuadPts[edgePos]) * vec2;
00640 Point rayref = innerQuadPts[edgePos] * vec1ref + (1.0
00641 - innerQuadPts[edgePos]) * vec2ref;
00642
00643 int nrIntersect = 0;
00644 Array<double> t;
00645 Point rayEnd = x + ray;
00646 curve.returnIntersect(x, rayEnd, nrIntersect, t);
00647
00648 double mu = 1.0;
00649 if (nrIntersect == 1)
00650 {
00651 mu = t[0];
00652 }
00653 else if (nrIntersect >= 2)
00654 {
00655
00656 integrals.resize(0);
00657 return;
00658 }
00659
00660
00661 Array<double> innerIntegrals(integrals.size());
00662 for (int i = 0; i < innerIntegrals.size(); i++)
00663 innerIntegrals[i] = 0.0;
00664
00665
00666 for (int rayPos = 0; rayPos < nrInnerQuadPts; rayPos++)
00667 {
00668 Point q = xref + mu * innerQuadPts[rayPos] * rayref;
00669
00670
00671 Array<double> funcVals;
00672 evaluateAllBasisFunctions(TriangleCell, q, funcVals);
00673
00674
00675 for (int i = 0; i < funcVals.size(); i++)
00676 innerIntegrals[i] += funcVals[i] * innerQuadWgts[rayPos]
00677 * innerQuadPts[rayPos];
00678 }
00679 for (int j = 0; j < integrals.size(); j++)
00680 integrals[j] += innerQuadWgts[edgePos] * mu * mu
00681 * innerIntegrals[j];
00682 }
00683
00684
00685 double det = ::fabs(vec1ref[0] * vec2ref[1] - vec2ref[0] * vec1ref[1]);
00686 for (int j = 0; j < integrals.size(); j++)
00687 integrals[j] *= det;
00688 } break;
00689
00690 case QuadCell:
00691 {
00692
00693 integrals.resize((order() + 1) * (order() + 1));
00694 for (int i = 0; i < integrals.size(); i++)
00695 integrals[i] = 0.0;
00696 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::integrateRegion: Initializing quadrature on a QuadCell")
00697
00698
00699 int nrInnerQuadPts = innerOrder + 3;
00700 nrInnerQuadPts = (nrInnerQuadPts + nrInnerQuadPts % 2) / 2;
00701 GaussLobatto1D innerQuad(nrInnerQuadPts, 0.0, 1.0);
00702
00703 Array<double> innerQuadPts(nrInnerQuadPts);
00704 Array<double> innerQuadWgts(nrInnerQuadPts);
00705 innerQuadPts = innerQuad.nodes();
00706 innerQuadWgts = innerQuad.weights();
00707
00708
00709 for (int xPos = 0; xPos < nrInnerQuadPts; xPos++)
00710 {
00711
00712 Point x_cur = x + innerQuadPts[xPos] * vec1;
00713 Point x_curref = xref + innerQuadPts[xPos] * vec1ref;
00714
00715 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::integrateRegion: Integrating at xPos " << xPos << " Ref point " << x_curref)
00716
00717
00718 int nrIntersect = 0;
00719 Array<double> t;
00720 Point x_curEnd = x_cur + vec2;
00721 curve.returnIntersect(x_cur, x_curEnd, nrIntersect, t);
00722 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::integrateRegion: Nr. of intersections along y: " << t.size())
00723
00724 double mu = 1.0;
00725 if (nrIntersect == 1)
00726 {
00727 mu = t[0];
00728 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::integrateRegion: Intersection parameter is " << mu)
00729 }
00730 else if (nrIntersect >= 2)
00731 {
00732
00733 integrals.resize(0);
00734 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::integrateRegion: More than one intersection -> Abort!")
00735
00736 return;
00737 }
00738
00739
00740 Array<double> innerIntegrals(integrals.size());
00741 for (int jj = 0; jj < innerIntegrals.size(); jj++)
00742 innerIntegrals[jj] = 0.0;
00743
00744 for (int yPos = 0; yPos < nrInnerQuadPts; yPos++)
00745 {
00746 Point y_curref = x_curref + mu * innerQuadPts[yPos] * vec2ref;
00747 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::integrateRegion: Integrating at yPos " << yPos << " Ref point " << y_curref)
00748
00749
00750 Array<double> funcVals;
00751 evaluateAllBasisFunctions(QuadCell, y_curref, funcVals);
00752
00753
00754 for (int jj = 0; jj < funcVals.size(); jj++)
00755 innerIntegrals[jj] += funcVals[jj] * innerQuadWgts[yPos] * mu;
00756 }
00757 for (int ii = 0; ii < integrals.size(); ii++)
00758 integrals[ii] += innerQuadWgts[xPos] * innerIntegrals[ii];
00759 }
00760
00761
00762 double det = ::fabs(vec1ref[0] * vec2ref[1] - vec2ref[0] * vec1ref[1]);
00763 for (int i = 0; i < integrals.size(); i++)
00764 integrals[i] *= det;
00765 SUNDANCE_MSG3(verb, tabs << "FeketeQuadrature::integrateRegion: Integration of region finished. Integral values " << integrals)
00766 } break;
00767
00768 default:
00769 #ifndef TRILINOS_7
00770 SUNDANCE_ERROR("FeketeQuadrature::integrateRegion() not available for cell type " << cellType)
00771 ;
00772 #else
00773 SUNDANCE_ERROR7("FeketeQuadrature::integrateRegion() not available for cell type " << cellType);
00774 #endif
00775 }
00776
00777 }
00778
00779 void FeketeQuadrature::evaluateAllBasisFunctions(const CellType cellType,
00780 const Point& q, Array<double>& result) const
00781 {
00782
00783 switch (cellType)
00784 {
00785 case TriangleCell:
00786 {
00787
00788
00789 if (!_hasBasisCoeffs)
00790 {
00791 FeketeTriangleQuadrature::computeBasisCoeffs(order(), _basisCoeffs);
00792 _hasBasisCoeffs = true;
00793 }
00794
00795
00796
00797
00798 int nFeketePts = sqrt((double) _basisCoeffs.size());
00799 result.resize(nFeketePts);
00800
00801
00802 Array<double> pkdVals(nFeketePts);
00803 FeketeTriangleQuadrature::evalPKDpolynomials(order(), q[0], q[1],
00804 &(pkdVals[0]));
00805
00806
00807
00808 char trans = 'N';
00809 double alpha = 1.;
00810 double beta = 0.;
00811 int inc = 1;
00812 ::dgemv_(&trans, &nFeketePts, &nFeketePts, &alpha, &(_basisCoeffs[0]),
00813 &nFeketePts, &(pkdVals[0]), &inc, &beta, &(result[0]), &inc);
00814 break;
00815 }
00816 case QuadCell:
00817 {
00818
00819 if (!_hasBasisCoeffs)
00820 {
00821 FeketeQuadQuadrature::computeBasisCoeffs(order(), _basisCoeffs);
00822 _hasBasisCoeffs = true;
00823 }
00824
00825
00826
00827
00828 int nFeketePts = sqrt((double) _basisCoeffs.size());
00829 result.resize(nFeketePts);
00830
00831
00832 Array<double> polVals(nFeketePts);
00833 FeketeQuadQuadrature::evalPolynomials(nFeketePts, q[0], q[1],
00834 &(polVals[0]));
00835
00836
00837
00838 char trans = 'N';
00839 double alpha = 1.;
00840 double beta = 0.;
00841 int inc = 1;
00842 ::dgemv_(&trans, &nFeketePts, &nFeketePts, &alpha, &(_basisCoeffs[0]),
00843 &nFeketePts, &(polVals[0]), &inc, &beta, &(result[0]), &inc);
00844 break;
00845 }
00846 default:
00847 #ifndef TRILINOS_7
00848 SUNDANCE_ERROR("FeketeQuadrature::evaluateAllBasisFunctions() not available for cell type " << cellType)
00849 ;
00850 #else
00851 SUNDANCE_ERROR7("FeketeQuadrature::evaluateAllBasisFunctions() not available for cell type " << cellType);
00852 #endif
00853 }
00854 }