306 lines
8.1 KiB
C++
306 lines
8.1 KiB
C++
#include "stdafx.h"
|
|
|
|
//***************************************************************************************************
|
|
// 坡屋顶求体积
|
|
//***************************************************************************************************
|
|
AcGePoint3d GetIntegralPointX(AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
|
|
{
|
|
AcGePoint3d ptRet;
|
|
//modified by fxj 比较浮点数大小用两个数做差与1E-6比
|
|
if(p1.x - p2.x > 1E-6 && p1.x - p3.x < -1E-6 || p1.x - p3.x > 1E-6 && p1.x - p2.x < -1E-6)
|
|
ptRet = p1;
|
|
else if(p2.x - p1.x > 1E-6 && p2.x - p3.x < -1E-6 || p2.x - p3.x > 1E-6 && p2.x - p1.x < -1E-6)
|
|
ptRet = p2;
|
|
else if(p3.x - p1.x > 1E-6 && p3.x - p2.x < -1E-6 || p3.x - p2.x > 1E-6 && p3.x - p1.x < -1E-6)
|
|
ptRet = p3;
|
|
else{
|
|
if(fabs(p1.x - p2.x) < 1E-6){
|
|
if(p1.y > p2.y)
|
|
ptRet = p1;
|
|
else
|
|
ptRet = p2;
|
|
}
|
|
else if(fabs(p1.x - p3.x) < 1E-6){
|
|
if(p1.y > p3.y)
|
|
ptRet = p1;
|
|
else
|
|
ptRet = p3;
|
|
}
|
|
else if(fabs(p2.x - p3.x) < 1E-6){
|
|
if(p2.y > p3.y)
|
|
ptRet = p2;
|
|
else
|
|
ptRet = p3;
|
|
}
|
|
}
|
|
return ptRet;
|
|
}
|
|
|
|
|
|
double GetIntegralX(AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
|
|
{
|
|
double dRet(0);
|
|
AcGePoint3d pt0, pt1, pt2;
|
|
pt0 = GetIntegralPointX(p1, p2, p3);
|
|
if(pt0 == p1){
|
|
if(p2.x < p3.x){
|
|
pt1 = p2;
|
|
pt2 = p3;
|
|
}
|
|
else{
|
|
pt1 = p3;
|
|
pt2 = p2;
|
|
}
|
|
}
|
|
else if(pt0 == p2){
|
|
if(p1.x < p3.x){
|
|
pt1 = p1;
|
|
pt2 = p3;
|
|
}
|
|
else{
|
|
pt1 = p3;
|
|
pt2 = p1;
|
|
}
|
|
}
|
|
else if(pt0 == p3){
|
|
if(p1.x < p2.x){
|
|
pt1 = p1;
|
|
pt2 = p2;
|
|
}
|
|
else{
|
|
pt1 = p2;
|
|
pt2 = p1;
|
|
}
|
|
}
|
|
|
|
double x0 = pt0.x;
|
|
double x1 = pt1.x;
|
|
double x2 = pt2.x;
|
|
double y0 = pt0.y;
|
|
double y1 = pt1.y;
|
|
double y2 = pt2.y;
|
|
|
|
double d1 = x0*x0-x1*x1;
|
|
double d2 = x2*x2-x0*x0;
|
|
double d3 = x0*x0*x0-x1*x1*x1;
|
|
double d4 = x2*x2*x2-x0*x0*x0;
|
|
|
|
|
|
double d6 = (y1-y2)/(x1-x2);
|
|
|
|
double nTmp3 = 1.0/3.0;
|
|
|
|
|
|
//modified by fxj 用三个点的顺逆时针判断倒正三角
|
|
XPoly2D poly(1);
|
|
poly.AppendNode(pt0);
|
|
poly.AppendNode(pt1);
|
|
poly.AppendNode(pt2);
|
|
double polyArea=poly.PathArea();
|
|
if(polyArea<-1E-6){//正三角形
|
|
double nTmp = 0;
|
|
if(fabs(x0-x1) > 1E-6 && fabs(x1 - x2) > 1E-6)//有2个点x坐标相同，只积分一次
|
|
{
|
|
double d5 = (y0-y1)/(x0-x1);
|
|
nTmp = 0.5*d1*(y0 - y1 - d5*x0 + d6*x1) + nTmp3*d3*(d5 - d6);
|
|
}
|
|
|
|
if(fabs(x1 - x2) > 1E-6 && fabs(x2 - x0) > 1E-6)
|
|
{
|
|
double d7 = (y2-y0)/(x2-x0);
|
|
nTmp += 0.5*d2*(y2 - y1 - d7*x2 + d6*x1) + nTmp3*d4*(d7 - d6);
|
|
}
|
|
dRet = nTmp;
|
|
}
|
|
else if(polyArea>1E-6){ //倒三角形
|
|
double nTmp = 0;
|
|
if(fabs(x0-x1) > 1E-6 && fabs(x1 - x2) > 1E-6)
|
|
{
|
|
double d5 = (y0-y1)/(x0-x1);
|
|
nTmp = 0.5*d1*(y1 - y0 - d6*x1 + d5*x0) + nTmp3*d3*(d6 - d5);
|
|
}
|
|
if(fabs(x1 - x2) > 1E-6 && fabs(x2 - x0) > 1E-6)
|
|
{
|
|
double d7 = (y2-y0)/(x2-x0);
|
|
nTmp += 0.5*d2*(y1 - y2 - d6*x1 + d7*x2) + nTmp3*d4*(d6 - d7);
|
|
}
|
|
dRet = nTmp;
|
|
}
|
|
return dRet;
|
|
}
|
|
|
|
|
|
|
|
double GetIntegralY( AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
|
|
{
|
|
AcGePoint3d ptTmp1, ptTmp2, ptTmp3;
|
|
|
|
ptTmp1.x = p1.y;
|
|
ptTmp1.y = p1.x;
|
|
ptTmp1.z = 0;
|
|
|
|
ptTmp2.x = p2.y;
|
|
ptTmp2.y = p2.x;
|
|
ptTmp2.z = 0;
|
|
ptTmp3.x = p3.y;
|
|
ptTmp3.y = p3.x;
|
|
ptTmp3.z = 0;
|
|
return GetIntegralX(ptTmp1, ptTmp2, ptTmp3);
|
|
|
|
}
|
|
|
|
|
|
//函数：通过积分，得到空间三角形投影到XY平面的体积
|
|
// 入口参数：
|
|
|
|
//p1, p2, p3 代表 投影三角形的3个顶点坐标
|
|
//Z0：积分水平基准面标高
|
|
//返回：体积积分
|
|
|
|
//modified by fxj
|
|
double TUTIL_IMPEXP GetTriangleVolumeByIntegral(AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3,double Z0=0)
|
|
{
|
|
p1.z-=Z0;
|
|
p2.z-=Z0;
|
|
p3.z-=Z0;
|
|
|
|
double dRet(0);
|
|
|
|
AcGeVector3d vt1(p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);
|
|
AcGeVector3d vt2(p3.x - p2.x, p3.y - p2.y, p3.z - p2.z);
|
|
AcGeVector3d vtNormal = vt1.crossProduct(vt2);
|
|
vtNormal.normalize();
|
|
AcGePlane plane(p1, vtNormal); //投影三角形的平面
|
|
double dA(0), dB(0), dC(0), dD(0); //dA, dB, dC,dD代表 aX + bY + cZ = D 空间三角形的平面方程
|
|
plane.getCoefficients(dA, dB, dC, dD);
|
|
dD = dD*(-1);
|
|
AcGeVector3d vtNormal2 = plane.normal();
|
|
|
|
BOOL bDzero=FALSE;
|
|
if(fabs(dD)>1E-6)//平面方程转化为 aX + bY + cZ = 1 形式
|
|
{
|
|
dA = dA / dD ;
|
|
dB = dB / dD ;
|
|
dC = dC / dD ;
|
|
|
|
bDzero=FALSE;
|
|
}
|
|
else
|
|
bDzero=TRUE;
|
|
|
|
if(fabs(dC)<1E-6)
|
|
return dRet;
|
|
|
|
XPoly2D poly;
|
|
poly.AppendNode(p1);
|
|
poly.AppendNode(p2);
|
|
poly.AppendNode(p3);
|
|
poly.AppendNode(p1);
|
|
poly.SetElevation(0);
|
|
|
|
double nProjectArea=poly.PathArea();
|
|
double nCalcTmp = 1/dC*nProjectArea;
|
|
double nX = dA / dC * GetIntegralX(p1, p2, p3);
|
|
double nY = dB / dC *GetIntegralY(p1, p2, p3);
|
|
|
|
if(nProjectArea<-1E-6)
|
|
{
|
|
nX=-nX;
|
|
nY=-nY;
|
|
}//积分是有方向的,与数学上积分的概念有所不同
|
|
|
|
|
|
dRet= - nX - nY;
|
|
if(bDzero!=TRUE)
|
|
dRet += nCalcTmp;
|
|
return dRet;
|
|
}
|
|
//对实心实体的体积计算 add by fxj
|
|
double VolumeEntity(AcDbObjectId id)
|
|
{
|
|
TMasterProperties prop;
|
|
AcDbVoidPtrArray entSet;
|
|
double dRet=0;
|
|
AcDbEntity *pEnt=NULL;
|
|
acdbOpenObject(pEnt,id,AcDb::kForRead);
|
|
if(pEnt!=NULL)
|
|
{
|
|
int nOldMode=DocGetDrawMode();
|
|
DocSetDrawMode(DRAW_3D);
|
|
|
|
AcDbVoidPtrArray TriangleSet;
|
|
InitModel(prop, pEnt, NULL, entSet);
|
|
DocSetDrawMode(nOldMode);
|
|
AcGeVector3d vtNormalZ(0,0,1);
|
|
|
|
int nLen=entSet.length();
|
|
for(int n=0; n<nLen; n++) {
|
|
AcDbEntity *pEnt=(AcDbEntity *)entSet[n];
|
|
if(pEnt->isKindOf(AcDbFace::desc())) {
|
|
AcDbFace *pFace = AcDbFace::cast(pEnt);
|
|
AcGePoint3d p1,p2,p3,p4;
|
|
pFace->getVertexAt(0, p1);
|
|
pFace->getVertexAt(1, p2);
|
|
pFace->getVertexAt(2, p3);
|
|
pFace->getVertexAt(3, p4);
|
|
AcGeVector3d vt1(p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);
|
|
AcGeVector3d vt2(p3.x - p2.x, p3.y - p2.y, p3.z - p2.z);
|
|
AcGeVector3d vtNormal = vt1.crossProduct(vt2);
|
|
vtNormal.normalize();
|
|
if(p1 == p4)
|
|
{
|
|
TriangleSet.append(pEnt->clone());
|
|
}
|
|
else
|
|
{
|
|
//注意保证逆时针方向
|
|
AcDbFace *pFaceNew1 = new AcDbFace();
|
|
pFaceNew1->setVertexAt(0, p1);
|
|
pFaceNew1->setVertexAt(1, p2);
|
|
pFaceNew1->setVertexAt(2, p3);
|
|
pFaceNew1->setVertexAt(3, p1);
|
|
TriangleSet.append(pFaceNew1);
|
|
AcDbFace *pFaceNew2 = new AcDbFace();
|
|
pFaceNew2->setVertexAt(0, p1);
|
|
pFaceNew2->setVertexAt(1, p3);
|
|
pFaceNew2->setVertexAt(2, p4);
|
|
pFaceNew2->setVertexAt(3, p1);
|
|
TriangleSet.append(pFaceNew2);
|
|
}
|
|
}
|
|
}
|
|
ReleaseEntSet(entSet);
|
|
nLen=TriangleSet.length();
|
|
|
|
for(n=0; n<nLen; n++) {
|
|
AcDbEntity *pEnt=(AcDbEntity *)TriangleSet[n];
|
|
if(pEnt->isKindOf(AcDbFace::desc())) {
|
|
AcDbFace *pFace = AcDbFace::cast(pEnt);
|
|
AcGePoint3d p1,p2,p3;
|
|
pFace->getVertexAt(0, p1);
|
|
pFace->getVertexAt(1, p2);
|
|
pFace->getVertexAt(2, p3);
|
|
double nV=0;
|
|
|
|
nV=GetTriangleVolumeByIntegral(p1, p2, p3);//三角形积分计算
|
|
|
|
dRet +=nV;
|
|
}
|
|
}
|
|
ReleaseEntSet(TriangleSet);
|
|
|
|
pEnt->close();
|
|
}
|
|
return -dRet;//与法向量有关
|
|
}
|
|
void VolumeRoofCmd()
|
|
{
|
|
Entity ent;
|
|
RBList rbMask(RTDXF0, _T("3DSOLID"), NULL);
|
|
if(SelectEntity(_T("\n请选择三维实体<退出>:"), ent, rbMask) != RTNORM)
|
|
return;
|
|
AcDbObjectId id = NULL;
|
|
acdbGetObjectId(id, ent);
|
|
double dV = VolumeEntity(id);
|
|
adsout << _T("\n三维实体的体积为：") << dV*1e-9 << _T("立方米\n");
|
|
} |