有限元框架中的数值积分部分。 #ifndef _SDL_MATHS_FEM_COMPONENT_INTEGRATION_HPP_ #define _SDL_MATHS_FEM_COMPONENT_INTEGRATION_HPP_ SDL_MATHS_FEM_BEGIN typedef Functor<FENumberType, TYPELIST_1(FENumberType)> FEFunction1DType; typedef Functor<FENumberType, TYPELIST_2(FENumberType, FENumberType)> FEFunction2DType; /* FENumberType TriangleIntegrate1(FEFunction2DType fun, const FEVertexListType * pvl) FENumberType TriangleIntegrate2(FEFunction2DType fun, const FEVertexListType * pvl) FENumberType TriangleIntegrate3(FEFunction2DType fun, const FEVertexListType * pvl) FENumberType TriangleIntegrate1(FEFunction2DType fun) FENumberType TriangleIntegrate2(FEFunction2DType fun) FENumberType TriangleIntegrate3(FEFunction2DType fun) FENumberType QuadrilateralIntegrate3(FEFunction2DType fun) FENumberType QuadrilateralIntegrate5(FEFunction2DType fun) FENumberType QuadrilateralIntegrate7(FEFunction2DType fun) #define QuadrilateralIntegrate1 QuadrilateralIntegrate3 #define QuadrilateralIntegrate2 QuadrilateralIntegrate3 #define QuadrilateralIntegrate4 QuadrilateralIntegrate5 #define QuadrilateralIntegrate6 QuadrilateralIntegrate7 FENumberType LineIntegrate3(FEFunction1DType fun, FENumberType from = -1.0, FENumberType to = 1.0) EENumberType LineIntegrate5(FEFunction1DType fun, FENumberType from = -1.0, FENumberType to = 1.0) FENumberType LineIntegrate7(FEFunction1DType fun, FENumberType from = -1.0, FENumberType to = 1.0) FENumberType LineIntegrate3(FEFunction2DType fun, FEFunction1DType y, FENumberType from = -1.0, FENumberType to = 1.0) FENumberType LineIntegrate5(FEFunction2DType fun, FEFunction1DType y, FENumberType from = -1.0, FENumberType to = 1.0) FENumberType LineIntegrate7(FEFunction2DType fun, FEFunction1DType y, FENumberType from = -1.0, FENumberType to = 1.0) #define LineIntegrate1 LineIntegrate3 #define LineIntegrate2 LineIntegrate3 #define LineIntegrate4 LineIntegrate5 #define LineIntegrate6 LineIntegrate7 */ inline FENumberType TriangleIntegrate1(FEFunction2DType fun, const FEVertexListType * pvl) { RT_CONDITION(pvl != NULL); RT_CONDITION(pvl->size() == 3); FENumberType area = CalculateArea<FENumberType>(*pvl); FEVertexType center = GetGravityCenter<FENumberType>(*pvl); return area * fun(center.x, center.y); } inline FENumberType TriangleIntegrate2(FEFunction2DType fun, const FEVertexListType * pvl) { RT_CONDITION(pvl != NULL); RT_CONDITION(pvl->size() == 3); FENumberType area = CalculateArea<FENumberType>(*pvl); FENumberType value = 0.0; for (Int32 i = 0; i < 3; ++i) { value += fun(((*pvl)[i].x+(*pvl)[(i+1)%3].x)/2, ((*pvl)[i].y+(*pvl)[(i+1)%3].y)/2); } return area / 3 * value; } inline FENumberType TriangleIntegrate3(FEFunction2DType fun, const FEVertexListType * pvl) { RT_CONDITION(pvl != NULL); RT_CONDITION(pvl->size() == 3); FENumberType area = CalculateArea<FENumberType>(*pvl); FEVertexType center = GetGravityCenter<FENumberType>(*pvl); FENumberType value = 27 * fun(center.x, center.y); for (Int32 i = 0; i < 3; ++i) { value += 3 * fun((*pvl)[i].x, (*pvl)[i].y); value += 8 * fun(((*pvl)[i].x+(*pvl)[(i+1)%3].x)/2, ((*pvl)[i].y+(*pvl)[(i+1)%3].y)/2); } return area / 60 * value; } inline FENumberType TriangleIntegrate1(FEFunction2DType fun) { static const FEVertexType vla[] = {FEVertexType(1,0), FEVertexType(0,1), FEVertexType(0, 0)}; static const FEVertexListType vl(vla, vla+3); return TriangleIntegrate1(fun, &vl); } inline FENumberType TriangleIntegrate2(FEFunction2DType fun) { static const FEVertexType vla[] = {FEVertexType(1,0), FEVertexType(0,1), FEVertexType(0, 0)}; static const FEVertexListType vl(vla, vla+3); return TriangleIntegrate2(fun, &vl); } inline FENumberType TriangleIntegrate3(FEFunction2DType fun) { static const FEVertexType vla[] = {FEVertexType(1,0), FEVertexType(0,1), FEVertexType(0, 0)}; static const FEVertexListType vl(vla, vla+3); return TriangleIntegrate3(fun, &vl); } inline FENumberType QuadrilateralIntegrate3(FEFunction2DType fun) { static const FENumberType gauss_point[] = {-0.5773502692, 0.5773502692}; static const FENumberType gauss_weight[] = {1.0, 1.0}; FENumberType result = 0.0; for (Int32 i = 0; i < 2; ++i) for (Int32 j = 0; j < 2; ++j) { result += gauss_weight[i] * gauss_weight[j] * fun(gauss_point[i], gauss_point[j]); } return result; } inline FENumberType QuadrilateralIntegrate5(FEFunction2DType fun) { static const FENumberType gauss_point[] = {-0.7745966692, 0, 0.7745966692}; static const FENumberType gauss_weight[] = {0.5555555556, 0.8888888889, 0.5555555556}; FENumberType result = 0.0; for (Int32 i = 0; i < 3; ++i) for (Int32 j = 0; j < 3; ++j) { result += gauss_weight[i] * gauss_weight[j] * fun(gauss_point[i], gauss_point[j]); } return result; } inline FENumberType QuadrilateralIntegrate7(FEFunction2DType fun) { static const FENumberType gauss_point[] = {-0.8611363116,-0.3399810436, 0.3399810436, 0.8611363116}; static const FENumberType gauss_weight[] = {0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451}; FENumberType result = 0.0; for (Int32 i = 0; i < 4; ++i) for (Int32 j = 0; j < 4; ++j) { result += gauss_weight[i] * gauss_weight[j] * fun(gauss_point[i], gauss_point[j]); } return result; } #define QuadrilateralIntegrate1 QuadrilateralIntegrate3 #define QuadrilateralIntegrate2 QuadrilateralIntegrate3 #define QuadrilateralIntegrate4 QuadrilateralIntegrate5 #define QuadrilateralIntegrate6 QuadrilateralIntegrate7 inline FENumberType LineIntegrate3(FEFunction1DType fun, FENumberType from = -1.0, FENumberType to = 1.0) { static const FENumberType gauss_point[] = {-0.5773502692, 0.5773502692}; static const FENumberType gauss_weight[] = {1.0, 1.0}; FENumberType result = 0.0; FENumberType d = (to - from)/2; FENumberType s = (from + to)/2; for (Int32 i = 0; i < 2; ++i) { result += gauss_weight[i] * fun(d*gauss_point[i]+s); } return d * result; } inline FENumberType LineIntegrate5(FEFunction1DType fun, FENumberType from = -1.0, FENumberType to = 1.0) { static const FENumberType gauss_point[] = {-0.7745966692, 0, 0.7745966692}; static const FENumberType gauss_weight[] = {0.5555555556, 0.8888888889, 0.5555555556}; FENumberType result = 0.0; FENumberType d = (to - from)/2; FENumberType s = (from + to)/2; for (Int32 i = 0; i < 3; ++i) { result += gauss_weight[i] * fun(d*gauss_point[i]+s); } return d * result; } inline FENumberType LineIntegrate7(FEFunction1DType fun, FENumberType from = -1.0, FENumberType to = 1.0) { static const FENumberType gauss_point[] = {-0.8611363116,-0.3399810436, 0.3399810436, 0.8611363116}; static const FENumberType gauss_weight[] = {0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451}; FENumberType result = 0.0; FENumberType d = (to - from)/2; FENumberType s = (from + to)/2; for (Int32 i = 0; i < 4; ++i) { result += gauss_weight[i] * fun(d*gauss_point[i]+s); } return d * result; } inline FENumberType LineIntegrate3(FEFunction2DType fun, FEFunction1DType y, FENumberType from = -1.0, FENumberType to = 1.0) { static const FENumberType gauss_point[] = {-0.5773502692, 0.5773502692}; static const FENumberType gauss_weight[] = {1.0, 1.0}; FENumberType result = 0.0; FENumberType d = (to - from)/2; FENumberType s = (from + to)/2; for (Int32 i = 0; i < 2; ++i) { FENumberType x = d * gauss_point[i] + s; result += gauss_weight[i] * fun(x, y(x)); } return d * result; } inline FENumberType LineIntegrate5(FEFunction2DType fun, FEFunction1DType y, FENumberType from = -1.0, FENumberType to = 1.0) { static const FENumberType gauss_point[] = {-0.7745966692, 0, 0.7745966692}; static const FENumberType gauss_weight[] = {0.5555555556, 0.8888888889, 0.5555555556}; FENumberType result = 0.0; FENumberType d = (to - from)/2; FENumberType s = (from + to)/2; for (Int32 i = 0; i < 3; ++i) { FENumberType x = d * gauss_point[i] + s; result += gauss_weight[i] * fun(x, y(x)); } return d * result; } inline FENumberType LineIntegrate7(FEFunction2DType fun, FEFunction1DType y, FENumberType from = -1.0, FENumberType to = 1.0) { static const FENumberType gauss_point[] = {-0.8611363116,-0.3399810436, 0.3399810436, 0.8611363116}; static const FENumberType gauss_weight[] = {0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451}; FENumberType result = 0.0; FENumberType d = (to - from)/2; FENumberType s = (from + to)/2; for (Int32 i = 0; i < 4; ++i) { FENumberType x = d * gauss_point[i] + s; result += gauss_weight[i] * fun(x, y(x)); } return d * result; } #define LineIntegrate1 LineIntegrate3 #define LineIntegrate2 LineIntegrate3 #define LineIntegrate4 LineIntegrate5 #define LineIntegrate6 LineIntegrate7 SDL_MATHS_FEM_END #endif