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 "SundanceClosedNewtonCotes.hpp"
00032 #include "SundanceTriangleQuadrature.hpp"
00033 #include "SundanceTetQuadrature.hpp"
00034
00035 using namespace Sundance;
00036 using namespace Teuchos;
00037
00038 ClosedNewtonCotes::ClosedNewtonCotes(int order)
00039 : QuadratureFamilyBase(order)
00040 {
00041 TEUCHOS_TEST_FOR_EXCEPTION(!(order==2 || order==3), std::runtime_error, "order "
00042 << order << " not supported by ClosedNewtonCotes");
00043 }
00044
00045 XMLObject ClosedNewtonCotes::toXML() const
00046 {
00047 XMLObject rtn("ClosedNewtonCotes");
00048 rtn.addAttribute("order", Teuchos::toString(order()));
00049 return rtn;
00050 }
00051
00052
00053 void ClosedNewtonCotes::getLineRule(Array<Point>& quadPoints,
00054 Array<double>& quadWeights) const
00055 {
00056 Array<double> x;
00057 if (order()==2)
00058 {
00059 x = tuple(0.0, 1.0);
00060 quadWeights = tuple(0.5, 0.5);
00061 }
00062 else
00063 {
00064 x = tuple(0.0, 0.5, 1.0);
00065 quadWeights = tuple(1.0/6.0, 2.0/3.0, 1.0/6.0);
00066 }
00067 quadPoints.resize(x.size());
00068
00069 for (int i=0; i<x.length(); i++)
00070 {
00071 quadPoints[i] = Point(x[i]);
00072 }
00073 }
00074
00075 void ClosedNewtonCotes::getTriangleRule(Array<Point>& quadPoints,
00076 Array<double>& quadWeights) const
00077 {
00078 if (order()==2)
00079 {
00080 quadPoints=tuple(Point(0.0, 0.0), Point(1.0, 0.0), Point(0.0, 1.0));
00081 quadWeights=tuple(1.0/3.0, 1.0/3.0, 1.0/3.0);
00082 }
00083 else
00084 {
00085 quadPoints = tuple(Point(0.0, 0.0), Point(0.5, 0.0), Point(1.0, 0.0),
00086 Point(0.5, 0.5), Point(0.0, 1.0), Point(0.0, 0.5),
00087 Point(1.0/3.0, 1.0/3.0));
00088 quadWeights = tuple(3.0/60.0, 8.0/60.0, 3.0/60.0,
00089 8.0/60.0, 3.0/60.0, 8.0/60.0,
00090 27.0/60.0);
00091 }
00092 for (int i=0; i<quadWeights.size(); i++)
00093 {
00094 quadWeights[i] *= 0.5;
00095 }
00096 }
00097