#include "stdafx.h"

#include "DbSlopeRoof.h"


//***************************************************************************************************
//                                        坡屋顶求体积
//***************************************************************************************************
AcGePoint3d GetIntegralPointX(AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
{
    AcGePoint3d ptRet;
    if(p1.x > p2.x && p1.x < p3.x || p1.x > p3.x && p1.x < p2.x)
        ptRet = p1;
    else if(p2.x > p1.x && p2.x < p3.x || p2.x > p3.x && p2.x < p1.x)
        ptRet = p2;
    else if(p3.x > p1.x && p3.x < p2.x || p3.x > p2.x && p3.x < p1.x)
        ptRet = p3;
    else{
        if(fabs(p1.x - p2.x) < 0.0001){
            if(p1.y > p2.y)
                ptRet = p1;
            else
                ptRet = p2;
        }
        else if(fabs(p1.x - p3.x) < 0.0001){
            if(p1.y > p3.y)
                ptRet = p1;
            else
                ptRet = p3;
        }
        else if(fabs(p2.x - p3.x) < 0.0001){
            if(p2.y > p3.y)
                ptRet = p2;
            else
                ptRet = p3;
        }
    }
    return ptRet;
}
    
double GetIntegralX(double dA, double dB, double dC, 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 d5 = (y0-y1)/(x0-x1);
    double d6 = (y1-y2)/(x1-x2);
    double d7 = (y2-y0)/(x2-x0);
    double nTmp3 = 1.0/3.0;

    if(pt0.y > pt1.y || pt0.y > pt2.y){//正三角形
        double nTmp = 0;
        if(fabs(x0-x1) > 0.00001 && fabs(x1 - x2) > 0.00001)//有2个点x坐标相同，只积分一次
            nTmp = 0.5*d1*(y0 - y1 - d5*x0 + d6*x1) + nTmp3*d3*(d5 - d6);
        if(fabs(x1 - x2) > 0.00001 && fabs(x2 - x0) > 0.00001)
            nTmp += 0.5*d2*(y2 - y1 - d7*x2 + d6*x1) + nTmp3*d4*(d7 - d6);
        dRet = nTmp;
    }
    else{                             //倒三角形
        double nTmp = 0;
        if(fabs(x0-x1) > 0.00001 && fabs(x1 - x2) > 0.00001)
            nTmp = 0.5*d1*(y1 - y0 - d6*x1 + d5*x0) + nTmp3*d3*(d6 - d5);
        if(fabs(x1 - x2) > 0.00001 && fabs(x2 - x0) > 0.00001)
            nTmp += 0.5*d2*(y1 - y2 - d6*x1 + d7*x2) + nTmp3*d4*(d6 - d7);
        dRet = nTmp;
    }
    return dRet;
}



double GetIntegralY(double dA, double dB, double dC, 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(dA, dB, dC, ptTmp1, ptTmp2, ptTmp3);

}


//函数：通过积分，得到空间三角形投影到XY平面的体积
//                入口参数：
//dA, dB, dC 代表 aX + bY + cZ = 1 空间三角形的平面方程
//p1, p2, p3 代表 投影三角形的3个顶点坐标
//返回：体积
double GetTriangleVolumeByIntegral(double dA, double dB, double dC, AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
{
    p1.x = p1.x / 1000; //转化为米
    p2.x = p2.x / 1000; //转化为米
    p3.x = p3.x / 1000; //转化为米
    p1.y = p1.y / 1000; //转化为米
    p2.y = p2.y / 1000; //转化为米
    p3.y = p3.y / 1000; //转化为米
    p1.z = p1.z / 1000; //转化为米
    p2.z = p2.z / 1000; //转化为米
    p3.z = p3.z / 1000; //转化为米
    double dRet(0);
    XPoly2D poly;
    poly.AppendNode(p1);
    poly.AppendNode(p2);
    poly.AppendNode(p3);
    poly.AppendNode(p1);
    double nProjectArea = poly.Area();
    double nCalcTmp = 1/dC*nProjectArea;
    double nX = dA / dC * GetIntegralX(dA, dB, dC, p1, p2, p3);
    double nY = dB / dC *GetIntegralY(dA, dB, dC, p1, p2, p3);

    dRet = nCalcTmp - nX - nY;
    return dRet;
}

double GetTriangleVolumeByIntegral(double dA, AcGePoint3d p1, AcGePoint3d p2, AcGePoint3d p3)
{
	double dHeight = ((p1.z - dA) + (p2.z - dA) + (p3.z - dA)) / 3;
	dHeight = dHeight / 1000;
	p1.x = p1.x / 1000; //转化为米
    p2.x = p2.x / 1000; //转化为米
    p3.x = p3.x / 1000; //转化为米
    p1.y = p1.y / 1000; //转化为米
    p2.y = p2.y / 1000; //转化为米
    p3.y = p3.y / 1000; //转化为米
    p1.z = 0.0;
    p2.z = 0.0;
    p3.z = 0.0;
	
    double dRet(0);
    XPoly2D poly;
    poly.AppendNode(p1);
    poly.AppendNode(p2);
    poly.AppendNode(p3);
    poly.AppendNode(p1);
    double nProjectArea = poly.Area();
	dRet = nProjectArea * dHeight;
	return dRet;
}

double VolumeRoof2(AcDbObjectId id)
{
    double dRet(0);
    OPENOBJ_BEGIN(id, AcDb::kForRead, TDbSlopeRoof, pRoof);
    if(pRoof){
        XPoly2D poly;
        poly = pRoof->GetOutlinePoly();
		double dElev = pRoof->GetElevation();
        AcGeVector3d vtNormalZ(0,0,1);
        AcGePlane    planeXY(*poly[0], vtNormalZ); //确定投影面
        AcDbVoidPtrArray entSet;
        AcDbVoidPtrArray TriangleSet;
        XExplodeStream dc(entSet);
        pRoof->ExplodeOrDraw3d(dc);
        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(vtNormal.z < 0.00001 )
                    continue;
                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 = GetTriangleVolumeByIntegral(dElev, p1, p2, p3);
                if(nV > 0)
                    dRet += nV;
            }
        }
        ReleaseEntSet(TriangleSet);
    }
    OPENOBJ_END();


	//人字
    OPENOBJ_BEGIN(id, AcDb::kForRead, TDbDualSlopeRoof, pRoof);
    if(pRoof){
		double dRoofEle(0);
		DocSetDrawMode(DRAW_3D);
		AcDbExtents ext;
		pRoof->getGeomExtents(ext);
		double dElev = pRoof->GetElevation();
		DocSetDrawMode(DRAW_2D);
		dRoofEle = pRoof->GetElevation();
		//	dRoofEle = ext.minPoint().z; // sr 2012-9-27 屋顶带厚度，应该直接读取GetElevation()

        XPoly2D poly;
        pRoof->GetPolyPath(poly);
        AcGeVector3d vtNormalZ(0,0,1);
        AcGePlane    planeXY(*poly[0], vtNormalZ); //确定投影面
        AcDbVoidPtrArray entSet;
        AcDbVoidPtrArray TriangleSet;
        XExplodeStream dc(entSet);
        pRoof->ExplodeOrDraw3d(dc);
        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(vtNormal.z < 0.00001 )
                    continue;
                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 = GetTriangleVolumeByIntegral(dElev, p1, p2, p3);
                if(nV > 0)
                    dRet += nV;
            }
        }
        ReleaseEntSet(TriangleSet);
    }
    OPENOBJ_END();    

	//攒尖屋顶
    OPENOBJ_BEGIN(id, AcDb::kForRead, TDbCuspRoof, pRoof);
    if(pRoof){
		double dRoofEle(0);
		DocSetDrawMode(DRAW_3D);
		AcDbExtents ext;
		pRoof->getGeomExtents(ext);
		DocSetDrawMode(DRAW_2D);
		dRoofEle = pRoof->GetElevation();
	//	dRoofEle = ext.minPoint().z; // sr 2012-9-27 屋顶带厚度，应该直接读取GetElevation()
        XPoly2D poly;
        pRoof->GetPolyPath(poly);
        AcGeVector3d vtNormalZ(0,0,1);
        AcGePlane    planeXY(*poly[0], vtNormalZ); //确定投影面
        AcDbVoidPtrArray entSet;
        AcDbVoidPtrArray TriangleSet;
        XExplodeStream dc(entSet);
		double dElev = pRoof->GetElevation();
        pRoof->ExplodeOrDraw3d(dc);
        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(vtNormal.z < 0.00001 )
                    continue;
                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();
        double dA(0), dB(0), dC(0), dD(0); 
        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 = GetTriangleVolumeByIntegral(dElev, p1, p2, p3);
                if(nV > 0)
                    dRet += nV;
            }
        }
        ReleaseEntSet(TriangleSet);
    }
    OPENOBJ_END();

	return dRet;
}

double VolumeRoof(AcDbObjectId id)
{
    double dRet(0);
    OPENOBJ_BEGIN(id, AcDb::kForRead, TDbSlopeRoof, pRoof);
    if(pRoof){
        XPoly2D poly;
        poly = pRoof->GetOutlinePoly();
        AcGeVector3d vtNormalZ(0,0,1);
        AcGePlane    planeXY(*poly[0], vtNormalZ); //确定投影面
        AcDbVoidPtrArray entSet;
        AcDbVoidPtrArray TriangleSet;
        XExplodeStream dc(entSet);
        pRoof->ExplodeOrDraw3d(dc);
        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(vtNormal.z < 0.00001 )
                    continue;
                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();
        double dA(0), dB(0), dC(0), dD(0); 
        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);
                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); //投影三角形的平面
                plane.getCoefficients(dA, dB, dC, dD);
                dD = dD*(-1);
                AcGeVector3d vtNormal2 = plane.normal();
                p1 = p1.project(planeXY, vtNormalZ);
                p2 = p2.project(planeXY, vtNormalZ);
                p3 = p3.project(planeXY, vtNormalZ);
                dA = dA / dD * 1000;
                dB = dB / dD * 1000;
                dC = dC / dD * 1000;
                double nang = vtNormalZ.angleTo(vtNormal);
                double nV = GetTriangleVolumeByIntegral(dA, dB, dC, p1, p2, p3);
                if(nV > 0)
                    dRet += nV;
                if(nV < 0){//积分算法有错误，在特点的▽▼△▲下，体积《0
					double nV2 = GetTriangleVolumeByIntegral(dA, dB, dC, p1, p2, p3);
                    AcGePoint3d p11,p22,p33;
                    pFace->getVertexAt(0, p11);
                    pFace->getVertexAt(1, p22);
                    pFace->getVertexAt(2, p33);
                    
                    AcDbFace *pFaceNew = new AcDbFace();
                    pFaceNew->setVertexAt(0, p11);
                    pFaceNew->setVertexAt(1, p22);
                    pFaceNew->setVertexAt(2, p33);
                    pFaceNew->setVertexAt(3, p11);
                    AddToModelSpace(pFaceNew, Adesk::kTrue);
                }
            }
        }
        ReleaseEntSet(TriangleSet);
        //因为计算完体积有bug，投影到0标高的体积。所以减去下面的部分。
        double dArea = poly.Area()/1000000;
        dRet -= dArea * pRoof->GetElevation()/1000.0;
    }
    OPENOBJ_END();


	//人字
    OPENOBJ_BEGIN(id, AcDb::kForRead, TDbDualSlopeRoof, pRoof);
    if(pRoof){
		double dRoofEle(0);
		DocSetDrawMode(DRAW_3D);
		AcDbExtents ext;
		pRoof->getGeomExtents(ext);
		DocSetDrawMode(DRAW_2D);
		dRoofEle = pRoof->GetElevation();
		//	dRoofEle = ext.minPoint().z; // sr 2012-9-27 屋顶带厚度，应该直接读取GetElevation()

        XPoly2D poly;
        pRoof->GetPolyPath(poly);
        AcGeVector3d vtNormalZ(0,0,1);
        AcGePlane    planeXY(*poly[0], vtNormalZ); //确定投影面
        AcDbVoidPtrArray entSet;
        AcDbVoidPtrArray TriangleSet;
        XExplodeStream dc(entSet);
        pRoof->ExplodeOrDraw3d(dc);
        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(vtNormal.z < 0.00001 )
                    continue;
                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();
        double dA(0), dB(0), dC(0), dD(0); 
        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);
                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); //投影三角形的平面
                plane.getCoefficients(dA, dB, dC, dD);
                dD = dD*(-1);
                AcGeVector3d vtNormal2 = plane.normal();
                p1 = p1.project(planeXY, vtNormalZ);
                p2 = p2.project(planeXY, vtNormalZ);
                p3 = p3.project(planeXY, vtNormalZ);
                dA = dA / dD * 1000;
                dB = dB / dD * 1000;
                dC = dC / dD * 1000;
                double nang = vtNormalZ.angleTo(vtNormal);
                double nV = GetTriangleVolumeByIntegral(dA, dB, dC, p1, p2, p3);
                if(nV > 0)
                    dRet += nV;
                if(nV < 0){//积分算法有错误，在特点的▽▼△▲下，体积《0
                    AcGePoint3d p11,p22,p33;
                    pFace->getVertexAt(0, p11);
                    pFace->getVertexAt(1, p22);
                    pFace->getVertexAt(2, p33);
                    
                    AcDbFace *pFaceNew = new AcDbFace();
                    pFaceNew->setVertexAt(0, p11);
                    pFaceNew->setVertexAt(1, p22);
                    pFaceNew->setVertexAt(2, p33);
                    pFaceNew->setVertexAt(3, p11);
                    AddToModelSpace(pFaceNew, Adesk::kTrue);
                }
            }
        }
        ReleaseEntSet(TriangleSet);
        //因为计算完体积有bug，投影到0标高的体积。所以减去下面的部分。
        double dArea = poly.Area()/1000000;
        dRet -= dArea * dRoofEle/1000.0;
    }
    OPENOBJ_END();    

	//攒尖屋顶
    OPENOBJ_BEGIN(id, AcDb::kForRead, TDbCuspRoof, pRoof);
    if(pRoof){
		double dRoofEle(0);
		DocSetDrawMode(DRAW_3D);
		AcDbExtents ext;
		pRoof->getGeomExtents(ext);
		DocSetDrawMode(DRAW_2D);
		dRoofEle = pRoof->GetElevation();
	//	dRoofEle = ext.minPoint().z; // sr 2012-9-27 屋顶带厚度，应该直接读取GetElevation()
        XPoly2D poly;
        pRoof->GetPolyPath(poly);
        AcGeVector3d vtNormalZ(0,0,1);
        AcGePlane    planeXY(*poly[0], vtNormalZ); //确定投影面
        AcDbVoidPtrArray entSet;
        AcDbVoidPtrArray TriangleSet;
        XExplodeStream dc(entSet);
        pRoof->ExplodeOrDraw3d(dc);
        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(vtNormal.z < 0.00001 )
                    continue;
                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();
        double dA(0), dB(0), dC(0), dD(0); 
        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);
                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); //投影三角形的平面
                plane.getCoefficients(dA, dB, dC, dD);
                dD = dD*(-1);
                AcGeVector3d vtNormal2 = plane.normal();
                p1 = p1.project(planeXY, vtNormalZ);
                p2 = p2.project(planeXY, vtNormalZ);
                p3 = p3.project(planeXY, vtNormalZ);
                dA = dA / dD * 1000;
                dB = dB / dD * 1000;
                dC = dC / dD * 1000;
                double nang = vtNormalZ.angleTo(vtNormal);
                double nV = GetTriangleVolumeByIntegral(dA, dB, dC, p1, p2, p3);
                if(nV > 0)
                    dRet += nV;
                if(nV < 0){//积分算法有错误，在特点的▽▼△▲下，体积《0
                    AcGePoint3d p11,p22,p33;
                    pFace->getVertexAt(0, p11);
                    pFace->getVertexAt(1, p22);
                    pFace->getVertexAt(2, p33);
                    
                    AcDbFace *pFaceNew = new AcDbFace();
                    pFaceNew->setVertexAt(0, p11);
                    pFaceNew->setVertexAt(1, p22);
                    pFaceNew->setVertexAt(2, p33);
                    pFaceNew->setVertexAt(3, p11);
                    AddToModelSpace(pFaceNew, Adesk::kTrue);
                }
            }
        }
        ReleaseEntSet(TriangleSet);
        //因为计算完体积有bug，投影到0标高的体积。所以减去下面的部分。
        double dArea = poly.Area()/1000000;
        dRet -= dArea * dRoofEle/1000.0;
    }
    OPENOBJ_END();

	return dRet;
}

void GetRoofAllPoly(AcDbObjectId id, DList<XPoly2D> &lstPoly)
{
    double dRet(0);
    OPENOBJ_BEGIN(id, AcDb::kForRead, TDbPolyPathEnt, pRoof);
    if(pRoof){
        AcDbVoidPtrArray entSet;
        XExplodeStream dc(entSet);
        pRoof->ExplodeOrDraw3d(dc);
        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(vtNormal.z < 0.00001)
                    continue;
                XPoly2D *poly = new XPoly2D();
                poly->AppendNode(p1);
                poly->AppendNode(p2);
                poly->AppendNode(p3);
                poly->AppendNode(p4);
                if(p1 != p4)
                    poly->AppendNode(p1);
                lstPoly.Append(poly);
            }
        }
        ReleaseEntSet(entSet);
    }
    OPENOBJ_END();
    return;
}


void VolumeRoofCmd2()
{
    Entity entRoof;
    RBList rbMask(RTDXF0, _T("ATS_SLOPEROOF,ATS_DUALROOF,ATS_CUSPROOF"),  NULL);
	if(SelectEntity(_T("\n请选择坡屋顶<退出>:"), entRoof, rbMask) != RTNORM)
        return;
    AcDbObjectId id = NULL;
    acdbGetObjectId(id, entRoof);
	double dV2 = VolumeRoof2(id);
    adsout << _T("\n新计算坡屋顶的体积为：") << dV2 << _T("立方米\n");
}

void VolumeRoofCmd()
{
    Entity entRoof;
    RBList rbMask(RTDXF0, _T("ATS_SLOPEROOF,ATS_DUALROOF,ATS_CUSPROOF"),  NULL);
	if(SelectEntity(_T("\n请选择坡屋顶<退出>:"), entRoof, rbMask) != RTNORM)
        return;
    AcDbObjectId id = NULL;
    acdbGetObjectId(id, entRoof);
    double dV = VolumeRoof(id);
    adsout << _T("\n坡屋顶的体积为：") << dV << _T("立方米\n");
    DList<XPoly2D> lstPoly;
    GetRoofAllPoly(id, lstPoly);
    int len = lstPoly.Length();
    double dArea(0);
    for(int i=0;i<len;i++){
        XPoly2D poly = *lstPoly[i];
       // poly.Make();
        AcGePoint3d p1,p2,p3;
        p1 = *poly[0];
        p2 = *poly[1];
        p3 = *poly[2];
        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();
        AcGeMatrix3d  mat;
        mat.setToWorldToPlane(vtNormal);
        poly.TransformBy(mat);
        poly.SetElevation(0);
        double dTmp = poly.Area()/1000000;
        dArea += dTmp;
    }
    adsout << _T("\n坡屋顶的面积为：") << dArea << _T("平方米\n");
}