出现无法读取内存情况怎么办？

1354081891 · 发表于 2018-4-2 11:44:07

马上注册，结交更多好友，享用更多功能^_^

您需要登录才可以下载或查看，没有账号？立即注册

x

这个程序做的是一个逾渗的程序，请问一下各位大佬，这种的情况改怎么找到错误出现的地方呀？
#include<iostream>
#include<ctime>
#include<cstdlib>
#include<fstream>
#include<cmath>
#include<iomanip>
#include<time.h>
#define DELFLAG -1//用来标记已经删除的粒子位置
using namespace std;
const double PI=3.14;
const int N=8;//非凸多边形的边数
const double L=50.0;//容器边长
ofstream outfile("F:\\全局坐标系下多边形的顶点坐标");
ofstream outfile1("F:\\Density-probility非周期边界条件");
ofstream outfile2("F:\\Density-probility周期边界条件");
ofstream outfile3("F:\\T-Density-probility非周期边界条件");
ofstream outfile4("F:\\T-Density-probility周期边界条件");
class point
{
public:
double x,y,z;
point()
:x(0),y(0),z(0)
{

}
};
class trianle//三角形类
{
public:
point *vertexs;//三角形顶点类
point *ProjectionVertexs;//三角形顶点的投影
trianle();
~trianle();
};
class polygon
{
public:
point polar_center;//极心坐标
double *polar_radius;//极径
double *polar_angle;//极角
point *PolygonVertexs;//多边形的顶点
trianle *DivideTriangles;//从多边形划分出的三角形
double a_fa,bei_ta,ga_ma;
point NoramalVector;//平面法向量;
polygon();
~polygon();
};
trianle::trianle()
{
vertexs=new point[3];
ProjectionVertexs=new point[3];//三角形顶点的投影
}
trianle::~trianle()
{
delete[] vertexs;
delete[] ProjectionVertexs;
}
polygon::polygon()
{
polar_center.x=0;polar_center.y=0;polar_center.z=0;
polar_angle=new double[N];
polar_radius=new double[N];
PolygonVertexs=new point[N];
DivideTriangles=new trianle[N];
a_fa=0.0;
bei_ta=0.0;
ga_ma=0.0;
NoramalVector.x=0.0;
NoramalVector.y=0.0;
NoramalVector.z=0.0;
}
polygon::~polygon()
{
delete[] polar_angle;
delete[] polar_radius;
delete[] PolygonVertexs;
delete[] DivideTriangles;
}
class PolygonsPercalation
{
private:
enum {num_rows=1};//表示粒子的密度组数
int N1[num_rows];//不同粒子密度下的粒子数
double density[num_rows];//每组粒子数对应的密度
double unperiodpercolationprobability[num_rows]; //每组density下对应非周期边界的percolation概率=T1/T;
double periodpercolationprobability[num_rows];//每组Vf下对应周期边界的percolation概率=T2/T;
polygon* A;//在给定的体积分数下的原始球的类，记为A类
polygon* B;//将A类中与容器边壁相交的粒子记录到B类中，B类是一个活动的类，以后粒子间相交，也记录到B类中
polygon* A1;//A1=A-B;即将A类中与B类相同的粒子删除，剩下A类中的粒子记为A1类，也为活动的类
polygon* AuxB;//从A类中找出与其中一个边壁相交的B类粒子的周期性补偿粒子，记为AuxB
public:
PolygonsPercalation();
~PolygonsPercalation();
void GetDate();//获取相关数据
void GeneratePolygon(polygon&);//生成非凸多边形
void GetTheLocatedCoordinate(polygon&);//计算局部坐标系下的坐标
bool delete_convex_polygon(polygon&);//剔除凸多边形
void RandomDispersion(int&);//将获取的N个多边形随机分布在容器中
void DividePolygonToTriangles(polygon&);//将多边形划分成三角形的函数
void extrapolate_the_polygon_to_3_D(polygon&);//将非凸多边形扩展至三维空间中,并确定平面的法向量

void PercolationSimulation();//逾渗的模拟
bool IntersectParticleToPlane(polygon&, double, double, double, double); //判断粒子与容器的面相交
void DetectPercolation_CurrplaneToParticles(int&, double, double, double, double, bool&, bool&); //考察当前的三个(左，上，前)面和粒子相交问题
void PeriodAux(polygon&, int&);
bool PercolationDectect(double, double, double, double, polygon*, polygon*, int&, int&, int&, bool&); //检测粒子间是否percolation,以及B类粒子和对应的面是否有相交
bool IntersectParticlesJudegmentInaccurately(polygon&, polygon&); //判断两个粒子相交粗糙判定
bool IntersectParticlesJudegementAccurately(polygon&, polygon&);//两不共面三角形相交精确判定
};
PolygonsPercalation::PolygonsPercalation()
{
A=new polygon[25000];
B=new polygon[25000];
A1=new polygon[25000];
AuxB=new polygon[10000];
}
PolygonsPercalation::~PolygonsPercalation()
{
delete [] A;
delete [] B;
delete [] A1;
delete [] AuxB;
}
void PolygonsPercalation::GetTheLocatedCoordinate(polygon& P)//计算局部坐标系下的顶点坐标
{
//***********************************************************************计算在局部坐标系下多边形的顶点坐标
   for(int t=0;t<N;t++)
   {
   P.PolygonVertexs[t].x= P.polar_radius[t]*cos( P.polar_angle[t]);
   P.PolygonVertexs[t].y= P.polar_radius[t]*sin( P.polar_angle[t]);
   P.PolygonVertexs[t].z=0;
   }
}
bool PolygonsPercalation::delete_convex_polygon(polygon& P)//排除凸多边形
{
double flag[N];//用于判定多边形凹凸性的参数
for(int i=0;i<N;i++)
{
int a=i+1,b=i+2;
if(i!=(N-2)&&i!=(N-1))
{
flag[i]=( P.PolygonVertexs[a].x- P.PolygonVertexs[i].x)*( P.PolygonVertexs[b].y- P.PolygonVertexs[a].y)-
( P.PolygonVertexs[b].x-P.PolygonVertexs[a].x)*(P.PolygonVertexs[a].y-P.PolygonVertexs[i].y);
}
if(i==(N-2))
{
flag[i]=(P.PolygonVertexs[a].x-P.PolygonVertexs[i].x)*(P.PolygonVertexs[0].y-P.PolygonVertexs[a].y)-
(P.PolygonVertexs[0].x-P.PolygonVertexs[a].x)*(P.PolygonVertexs[a].y-P.PolygonVertexs[i].y);
}
if(i==(N-1))
{
flag[i]=(P.PolygonVertexs[0].x-P.PolygonVertexs[i].x)*(P.PolygonVertexs[1].y-P.PolygonVertexs[0].y)-
(P.PolygonVertexs[1].x-P.PolygonVertexs[0].x)*(P.PolygonVertexs[0].y-P.PolygonVertexs[i].y);
}
}
if(flag[0]<0||flag[1]<0||flag[2]<0||flag[3]<0||flag[3]<0||flag[4]<0||flag[5]<0||flag[6]<0||flag[7]<0)
return (1);//多边形为非凸时返回为1
else
return 0;//为凸边形时返回为0
}
void PolygonsPercalation::GeneratePolygon(polygon& P)
{
double min_L,max_L;//容器边的最小，最大值
      double max_yita=1,min_yita=0,min_A=0.5,average_A=1;//min_A,average_A是粒子的最小和平均粒径
      double max_yita1=1,min_yita1=0;//生成极角和极径的参数yita
      double temp;
      min_L=-0.5*L;max_L=0.5*L;
      P.polar_center.x=((((float)rand()/RAND_MAX))*(max_L-min_L)+min_L);//随机生成在x轴的极心坐标
      P.polar_center.y=((((float)rand()/RAND_MAX))*(max_L-min_L)+min_L);//随机生成在y轴的极心坐标
      P.polar_center.z=((((float)rand()/RAND_MAX))*(max_L-min_L)+min_L);//随机生成在z轴的极心坐标
      for(int i=0;i<N;i++)//生成极心和极角
      {
      P.polar_angle[i]=((((float)rand()/RAND_MAX))*(max_yita-min_yita)+min_yita)*(2*PI);
      P.polar_radius[i]=average_A+(2*((((float)rand()/RAND_MAX))*(max_yita-min_yita)+min_yita)-1)+min_A;
      }
      for(int i=0;i<N;i++)//将随机生成的极角按照从小到大的顺序排列
      {
      for(int j=0;j<N-i-1;j++)
      {
      if(P.polar_angle[j]>P.polar_angle[j+1])
      {temp=P.polar_angle[j];P.polar_angle[j]=P.polar_angle[j+1];P.polar_angle[j+1]=temp;}
      }
      }
//**********************************************************************随机地生成三个欧拉角
      double max_a_fa=PI,min_a_fa=-PI;
      double max_ga_ma=PI,min_ga_ma=-PI;
      double max_bei_ta=PI,min_bei_ta=0;

      P.a_fa=((((float)rand()/RAND_MAX))*(max_a_fa-min_a_fa)+min_a_fa);
      P.ga_ma=((((float)rand()/RAND_MAX))*(max_ga_ma-min_ga_ma)+min_ga_ma);
      P.bei_ta=((((float)rand()/RAND_MAX))*(max_bei_ta-min_bei_ta)+min_bei_ta);
}
void PolygonsPercalation::DividePolygonToTriangles(polygon& P)
{
for(int i=0;i<N;i++)
{
int a=i+1;
if(i!=(N-1))
{
P.DivideTriangles[i].vertexs[0].x=P.PolygonVertexs[i].x;
P.DivideTriangles[i].vertexs[0].y=P.PolygonVertexs[i].y;
P.DivideTriangles[i].vertexs[0].z=P.PolygonVertexs[i].z;

P.DivideTriangles[i].vertexs[1].x=P.PolygonVertexs[a].x;
P.DivideTriangles[i].vertexs[1].y=P.PolygonVertexs[a].y;
P.DivideTriangles[i].vertexs[1].z=P.PolygonVertexs[a].z;

P.DivideTriangles[i].vertexs[2].x=P.polar_center.x;
P.DivideTriangles[i].vertexs[2].y=P.polar_center.y;
P.DivideTriangles[i].vertexs[2].z=P.polar_center.z;
//向XOY面投影的三角形顶点坐标
P.DivideTriangles[i].ProjectionVertexs[0].x=P.PolygonVertexs[i].x;
P.DivideTriangles[i].ProjectionVertexs[0].y=P.PolygonVertexs[i].y;
P.DivideTriangles[i].ProjectionVertexs[0].z=0;

P.DivideTriangles[i].ProjectionVertexs[1].x=P.PolygonVertexs[a].x;
P.DivideTriangles[i].ProjectionVertexs[1].y=P.PolygonVertexs[a].y;
P.DivideTriangles[i].ProjectionVertexs[1].z=0;

P.DivideTriangles[i].ProjectionVertexs[2].x=P.polar_center.x;
P.DivideTriangles[i].ProjectionVertexs[2].y=P.polar_center.y;
P.DivideTriangles[i].ProjectionVertexs[2].z=0;
}
if(i==(N-1))
{
P.DivideTriangles[i].vertexs[0].x=P.PolygonVertexs[i].x;
P.DivideTriangles[i].vertexs[0].y=P.PolygonVertexs[i].y;
P.DivideTriangles[i].vertexs[0].z=P.PolygonVertexs[i].z;

P.DivideTriangles[i].vertexs[1].x=P.PolygonVertexs[0].x;
P.DivideTriangles[i].vertexs[1].y=P.PolygonVertexs[0].y;
P.DivideTriangles[i].vertexs[1].z=P.PolygonVertexs[0].z;

P.DivideTriangles[i].vertexs[2].x=P.polar_center.x;
P.DivideTriangles[i].vertexs[2].y=P.polar_center.y;
P.DivideTriangles[i].vertexs[2].z=P.polar_center.z;

//向XOY面投影的三角形顶点坐标
P.DivideTriangles[i].ProjectionVertexs[0].x=P.PolygonVertexs[i].x;
P.DivideTriangles[i].ProjectionVertexs[0].y=P.PolygonVertexs[i].y;
P.DivideTriangles[i].ProjectionVertexs[0].z=0;

P.DivideTriangles[i].ProjectionVertexs[1].x=P.PolygonVertexs[0].x;
P.DivideTriangles[i].ProjectionVertexs[1].y=P.PolygonVertexs[0].y;
P.DivideTriangles[i].ProjectionVertexs[1].z=0;

P.DivideTriangles[i].ProjectionVertexs[2].x=P.polar_center.x;
P.DivideTriangles[i].ProjectionVertexs[2].y=P.polar_center.y;
P.DivideTriangles[i].ProjectionVertexs[2].z=0;
}
}
}
void PolygonsPercalation::extrapolate_the_polygon_to_3_D(polygon& P)//将局部坐标系下的顶点坐标外推到三维空间中
{
//***********************************************************************计算在全局坐标系下的坐标，并将其输出到文件中"
for(int k=0;k<N;k++)
{
P.PolygonVertexs[k].x=(cos(P.ga_ma)*cos(P.bei_ta))*P.PolygonVertexs[k].x+((-sin(P.ga_ma)*cos(P.a_fa))-cos(P.ga_ma)*sin(P.bei_ta)*sin(P.a_fa))*P.PolygonVertexs[k].y
+(sin(P.ga_ma)*sin(P.a_fa)-cos(P.ga_ma)*sin(P.bei_ta)*cos(P.a_fa))*P.PolygonVertexs[k].z+P.polar_center.x;
//**********************************************************************************计算在全局坐标系下的x坐标
P.PolygonVertexs[k].y=(sin(P.ga_ma)*cos(P.bei_ta))*P.PolygonVertexs[k].x+(cos(P.ga_ma)*cos(P.ga_ma)-sin(P.ga_ma)*sin(P.bei_ta)*sin(P.a_fa))*P.PolygonVertexs[k].y+
(-cos(P.ga_ma)*sin(P.a_fa)-sin(P.ga_ma)*sin(P.bei_ta)*cos(P.a_fa))*P.PolygonVertexs[k].z+P.polar_center.y;
//**********************************************************************************计算在全局坐标系下的y坐标
P.PolygonVertexs[k].z=sin(P.bei_ta)*P.PolygonVertexs[k].x+(cos(P.bei_ta)*sin(P.a_fa))*P.PolygonVertexs[k].y+
(cos(P.bei_ta)*cos(P.a_fa))*P.PolygonVertexs[k].z+P.polar_center.z;
//**********************************************************************************计算在全局坐标系下的z坐标
}
/*
for(int k=0;k<N;k++)
{
int a=k+1;
if(a==N)
{
outfile<<"polygon{3,<"<<PolygonVertexs[k].x<<","<<PolygonVertexs[k].y<<","<<PolygonVertexs[k].z<<">, <"
      <<polar_center.x<<","<<polar_center.y<<","<<polar_center.z<<">, <"
      <<PolygonVertexs[0].x<<","<<PolygonVertexs[0].y<<","<<PolygonVertexs[0].z<<">}"<<endl;
}
else if(a!=N)
{
outfile<<"polygon{3,<"<<PolygonVertexs[k].x<<","<<PolygonVertexs[k].y<<","<<PolygonVertexs[k].z<<">, <"
            <<polar_center.x<<","<<polar_center.y<<","<<polar_center.z<<">, <"
      <<PolygonVertexs[a].x<<","<<PolygonVertexs[a].y<<","<<PolygonVertexs[a].z<<">}"<<endl;
}
}*/
double mo_liang;
//*********************************************************************获取平面的法向量
P.NoramalVector.x=(P.PolygonVertexs[0].y-P.PolygonVertexs[1].y)*(P.PolygonVertexs[2].z-P.PolygonVertexs[1].z)-
                     (P.PolygonVertexs[0].z-P.PolygonVertexs[1].z)*(P.PolygonVertexs[2].y-P.PolygonVertexs[1].y);
P.NoramalVector.y=(P.PolygonVertexs[0].z-P.PolygonVertexs[1].z)*(P.PolygonVertexs[2].x-P.PolygonVertexs[1].x)-
                     (P.PolygonVertexs[0].x-P.PolygonVertexs[1].x)*(P.PolygonVertexs[2].z-P.PolygonVertexs[1].z);
P.NoramalVector.z=(P.PolygonVertexs[0].x-P.PolygonVertexs[1].x)*(P.PolygonVertexs[2].y-P.PolygonVertexs[1].y)-
                     (P.PolygonVertexs[0].y-P.PolygonVertexs[1].y)*(P.PolygonVertexs[2].x-P.PolygonVertexs[1].x);
mo_liang=sqrt(pow(P.NoramalVector.x,2)+pow(P.NoramalVector.y,2)+pow(P.NoramalVector.z,2));
P.NoramalVector.x/=mo_liang;
P.NoramalVector.y/=mo_liang;
P.NoramalVector.z/=mo_liang;

}
void PolygonsPercalation::GetDate()
{
cout<<"请输入15组不同的容器中颗粒密度"<<endl;
for(int i=0;i<num_rows;i++)
{
cin>>density[i];//输入容器中颗粒的密度
}
}
void PolygonsPercalation::PercolationSimulation()
{
for(int i=0;i<num_rows;i++)
{
double T=0.0; //程序循环的次数，主要为了计算percolation probability
double T1=0.0; //记录非周期边界条件下的T次循环中达到percolation的循环次数
double T2=0.0; //记录周期边界条件下的T次循环中达到percolation的循环次数
N1[i]=int(density[i]*pow(L,3));//int()强制类型转换
while(T<200)
{
bool valid_1_rightplane=false; //bool型变量，申明粒子与非周期右边壁不相交
bool valid_2_rightplane=false; //bool型变量，申明粒子与周期右边壁不相交

bool valid_1_upplane=false;//bool型变量，申明粒子与非周期上边壁不相交
bool valid_2_upplane=false;//bool型变量，申明粒子与周期上边壁不相交

bool valid_1_frontplane=false;//bool型变量，申明粒子与非周期前边壁不相交
bool valid_2_frontplane=false;//bool型变量，申明粒子与周期前边壁不相交

RandomDispersion(N1[i]); //先在容器中随机分布粒子，粒子间允许重叠发生
//考察容器的右边面是否有粒子与其构成percolation
DetectPercolation_CurrplaneToParticles(N1[i],1.0, 0.0, 0.0, -0.5*L, valid_1_rightplane, valid_2_rightplane);
if (valid_2_rightplane==false)//表明右边面无效，则需要考察另外一个面，即上边面
{
//考察容器的上边面是否有粒子与其构成percolation
cout<<"Enter testing UpPlane!!!!!!"<<endl;
DetectPercolation_CurrplaneToParticles(N1[i], 0.0, 0.0, 1.0,-0.5*L, valid_1_upplane, valid_2_upplane);
}
if(valid_2_rightplane==false && valid_2_upplane==false)
{
//考察容器的前边面是否有粒子与其构成percolation
cout<<"Enter testing FrontPlane!!!!!!"<<endl;
DetectPercolation_CurrplaneToParticles(N1[i],0.0, 1.0, 0.0,-0.5*L, valid_1_frontplane, valid_2_frontplane);
}
if (valid_1_rightplane==true || valid_1_upplane==true || valid_1_frontplane==true)
{
cout<<"Successful!!!!!!!!"<<endl;
T1++;
}
if (valid_2_rightplane==true || valid_2_upplane==true || valid_2_frontplane==true)
{

cout<<"periodic boundary condtions Successful!!!!!!!!"<<endl;
T2++;
}
cout<<T<<endl;
T++;
if (((int)T%10)==0)//这个条件语句是为了得到模拟次数与渗流概率的对应关系
{
unperiodpercolationprobability[i]=T1/T; //计算非周期边界条件下的percolation概率
periodpercolationprobability[i]=T2/T; //计算周期边界条件下的percolation概率
outfile3<<"T="<<T<<endl<<"density="<<density[i]<<" "<<unperiodpercolationprobability[i]<<endl;
outfile4<<"T="<<T<<endl<<"density="<<density[i]<<" "<<periodpercolationprobability[i]<<endl;
}
}
unperiodpercolationprobability[i]=T1/T; //计算非周期边界条件下的percolation概率
periodpercolationprobability[i]=T2/T; //计算周期边界条件下的percolation概率
outfile1<<density[i]<<" "<<unperiodpercolationprobability[i]<<endl;
outfile2<<density[i]<<" "<<periodpercolationprobability[i]<<endl;
}
}
void PolygonsPercalation::RandomDispersion(int& N1)
{
int T;
for(int k=0;k<N1;k++)
{
do
{
GeneratePolygon(A[k]);
GetTheLocatedCoordinate(A[k]);
}while(delete_convex_polygon(A[k])==0);
extrapolate_the_polygon_to_3_D(A[k]);
      DividePolygonToTriangles(A[k]);
}
}
void PolygonsPercalation::DetectPercolation_CurrplaneToParticles (int& N2, double a0, double b0, double c0, double d0, bool& valplane1, bool& valplane2)
{
int currb_count=0; //记录添加到B类中的粒子
int a1_count=0; //记录添加到A1类中的粒子：A1=A-B
bool currplane=false; //考察当前的面是否有粒子与其相交
bool parperco=false; //bool型变量，申明当前的面与粒子间不发生渗流
int auxb_count=0; //记录添加到AuxB类中的粒子数量
for(int i=0;i<N2;i++)//然后判断A类中粒子是否与容器边壁相交，若相交，记录相交的粒子作为B类
{
if(IntersectParticleToPlane(A[i],a0,b0,c0,d0)==true)
{
currplane=true; //表面有粒子与当前的面相交，那么该面可以继续考察
B[currb_count]=A[i];
PeriodAux(B[currb_count], auxb_count);
currb_count++;
}
else
{
A1[a1_count]=A[i]; //将与边壁不相交的粒子记为A1类
a1_count++;
}
}
if (currplane==true)//表示currB类中存在粒子，即A类中存在粒子与当前的面相交
{
parperco=PercolationDectect(a0, b0, c0, d0+L, B, A1, currb_count, a1_count, auxb_count, valplane1); //判断B类中是否有粒子与当前面的对应面相交，若有则percolation成立,parperco=true；反之，parperco=false.
if (parperco==true)//表示当前的面与粒子间发生percolation,该考察面成立
{
valplane2=true; //记录该面有效
}
else
valplane2=false; //表明当前的考察面无法达到percolation,即无效
}
else
valplane2=false;//表明A类中没有粒子与当前的考察面相交，即无法达到percolation
}
bool PolygonsPercalation::IntersectParticleToPlane(polygon& oc, double a, double b, double c, double d)//判断粒子与面是否相交
{
double DistanceWithDirection[3];//有方向的距离
if(b==1.0&&a==0&&c==0)
{
for(int i=0;i<N;i++)
      {
DistanceWithDirection[0]=oc.DivideTriangles[i].vertexs[0].x*a+(oc.DivideTriangles[i].vertexs[0].y-d)*b+oc.DivideTriangles[i].vertexs[0].z*c;
      DistanceWithDirection[1]=oc.DivideTriangles[i].vertexs[1].x*a+(oc.DivideTriangles[i].vertexs[1].y-d)*b+oc.DivideTriangles[i].vertexs[1].z*c;
      DistanceWithDirection[2]=oc.DivideTriangles[i].vertexs[2].x*a+(oc.DivideTriangles[i].vertexs[2].y-d)*b+oc.DivideTriangles[i].vertexs[2].z*c;
      if(DistanceWithDirection[0]>0&&DistanceWithDirection[1]>0&&DistanceWithDirection[2]>0)
continue;//三角形位于坐标面的同侧，不相交继续执行
      if(DistanceWithDirection[0]<0&&DistanceWithDirection[1]<0&&DistanceWithDirection[2]<0)
continue;//三角形位于坐标面的同侧，不相交继续执行
      else
return true;//其他情况，即DiStanceWIthDirection三个正负不同，或者全为0，返回为真
      }
return false;//所有的多边形与坐标面都不相交返回为假
}
if(a==1&&b==0&&c==0)
{
for(int i=0;i<N;i++)
      {
DistanceWithDirection[0]=(oc.DivideTriangles[i].vertexs[0].x-d)*a+(oc.DivideTriangles[i].vertexs[0].y)*b+oc.DivideTriangles[i].vertexs[0].z*c;
      DistanceWithDirection[1]=(oc.DivideTriangles[i].vertexs[1].x-d)*a+(oc.DivideTriangles[i].vertexs[1].y)*b+oc.DivideTriangles[i].vertexs[1].z*c;
      DistanceWithDirection[2]=(oc.DivideTriangles[i].vertexs[2].x-d)*a+(oc.DivideTriangles[i].vertexs[2].y)*b+oc.DivideTriangles[i].vertexs[2].z*c;
      if(DistanceWithDirection[0]>0&&DistanceWithDirection[1]>0&&DistanceWithDirection[2]>0)
continue;
      if(DistanceWithDirection[0]<0&&DistanceWithDirection[1]<0&&DistanceWithDirection[2]<0)
continue;
      else
return true;
      }
return false;
}
if(c==1&&a==0&&b==0)
{
for(int i=0;i<N;i++)
      {
DistanceWithDirection[0]=oc.DivideTriangles[i].vertexs[0].x*a+(oc.DivideTriangles[i].vertexs[0].y)*b+(oc.DivideTriangles[i].vertexs[0].z-d)*c;
      DistanceWithDirection[1]=oc.DivideTriangles[i].vertexs[1].x*a+(oc.DivideTriangles[i].vertexs[1].y)*b+(oc.DivideTriangles[i].vertexs[1].z-d)*c;
      DistanceWithDirection[2]=oc.DivideTriangles[i].vertexs[2].x*a+(oc.DivideTriangles[i].vertexs[2].y)*b+(oc.DivideTriangles[i].vertexs[2].z-d)*c;
      if(DistanceWithDirection[0]>0&&DistanceWithDirection[1]>0&&DistanceWithDirection[2]>0)
         continue;
      if(DistanceWithDirection[0]<0&&DistanceWithDirection[1]<0&&DistanceWithDirection[2]<0)
         continue;
      else
return true;
      }
return false;
}
}
void PolygonsPercalation::PeriodAux(polygon& oc,int& aux_count)
{
int curraux_count=0; //记录当前考察的粒子需要周期性补偿几个粒子数量
point* curraux; //记录需要周期性补偿的粒子
curraux=new point[7];
bool I1=1, I2=1, I3=1, I4=1, I5=1, I6=1; //粒子与立方体6个面(1,2,3,4,5,6)相交的变量
I1=IntersectParticleToPlane(oc, 0.0, 1.0, 0.0, 0.5*L);
I2=IntersectParticleToPlane(oc, 0.0, 1.0, 0.0, -0.5*L);
I3=IntersectParticleToPlane(oc, 1.0, 0.0, 0.0, -0.5*L);
I4=IntersectParticleToPlane(oc, 1.0, 0.0, 0.0, 0.5*L);
I5=IntersectParticleToPlane(oc, 0.0, 0.0, 1.0, -0.5*L);
I6=IntersectParticleToPlane(oc, 0.0, 0.0, 1.0, 0.5*L);
if(I2==1 && I3==0 && I4==0 && I5==0 && I6==0)//如果粒子只与面2相交，补偿1个附加球
{
curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y-L;
curraux[0].z=oc.polar_center.z;
curraux_count=1;
}
else if(I3==1 && I1==0 && I2==0 && I5==0 && I6==0)//如果粒子只与面3相交，补偿1个附加球
{
curraux[0].x=oc.polar_center.x-L;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z;
curraux_count=1;
}
else if(I5==1 && I1==0 && I2==0 && I3==0 && I4==0)//如果粒子只与面5相交，补偿1个附加球
{
curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z-L;
curraux_count=1;
}
else if(I1==1 && I3==1 && I5==0 && I6==0)//如果粒子只与面2和面3相交A1E1，补偿3个附加球
{
curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y+L;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x-L;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x-L;
curraux[2].y=oc.polar_center.y+L;
curraux[2].z=oc.polar_center.z;
curraux_count=3;
}
else if(I2==1 && I3==1 && I5==0 && I6==0)//如果粒子只与面1和面3相交B1F1，补偿3个附加球
{
curraux[0].x=oc.polar_center.x-L;
curraux[0].y=oc.polar_center.y-L;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x-L;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y-L;
curraux[2].z=oc.polar_center.z;
curraux_count=3;
}
else if(I2==1 && I4==1 && I5==0 && I6==0)//如果粒子只与面2和面4相交C1G1，补偿3个附加球
{
curraux[0].x=oc.polar_center.x+L;
curraux[0].y=oc.polar_center.y-L;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x+L;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y-L;
curraux[2].z=oc.polar_center.z;
curraux_count=3;
}
else if(I1==1 && I5==1 && I3==0 && I4==0)//如果粒子只与面1和面5相交A1D1，补偿3个附加球
{
curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y+L;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z-L;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y+L;
curraux[2].z=oc.polar_center.z-L;
curraux_count=3;
}
else if(I3==1 && I5==1 && I1==0 && I2==0)//如果粒子只与面3和面5相交A1B1，补偿3个附加球
{

curraux[0].x=oc.polar_center.x-L;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z-L;

curraux[2].x=oc.polar_center.x-L;
curraux[2].y=oc.polar_center.y;
curraux[2].z=oc.polar_center.z-L;
curraux_count=3;
}
else if(I2==1 && I5==1 && I3==0 && I4==0)//如果粒子只与面2和面5相交B1C1，补偿3个附加球
{
curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y-L;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z-L;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y-L;
curraux[2].z=oc.polar_center.z-L;
curraux_count=3;
}
else if(I4==1 && I5==1 && I1==0 && I2==0)//如果粒子只与面4和面5相交C1D1，补偿3个附加球
{

curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z-L;

curraux[1].x=oc.polar_center.x+L;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x+L;
curraux[2].y=oc.polar_center.y;
curraux[2].z=oc.polar_center.z-L;
curraux_count=3;
}
else if(I3==1 && I6==1 && I1==0 && I2==0)//如果粒子只与面3和面6相交E1F1，补偿3个附加球
{

curraux[0].x=oc.polar_center.x-L;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z+L;

curraux[2].x=oc.polar_center.x-L;
curraux[2].y=oc.polar_center.y;
curraux[2].z=oc.polar_center.z+L;
curraux_count=3;
}
else if(I2==1 && I6==1 && I3==0 && I4==0) //如果粒子只与面2和面6相交F1G1，补偿3个附加球
{

curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z+L;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y-L;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y-L;
curraux[2].z=oc.polar_center.z+L;
curraux_count=3;
}
else if(I1==1 && I3==1 && I5==1)//如果粒子与面1，面3和面5相交类似于角A，补偿7个附加球
{

curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z-L;

curraux[1].x=oc.polar_center.x-L;
curraux[1].y=oc.polar_center.y;
curraux[1].y=oc.polar_center.z;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y+L;
curraux[2].z=oc.polar_center.z;

curraux[3].x=oc.polar_center.x-L;
curraux[3].y=oc.polar_center.y+L;
curraux[3].z=oc.polar_center.z;

curraux[4].x=oc.polar_center.x-L;
curraux[4].y=oc.polar_center.y;
curraux[4].z=oc.polar_center.z-L;

curraux[5].x=oc.polar_center.x;
curraux[5].y=oc.polar_center.y+L;
curraux[5].z=oc.polar_center.z-L;

curraux[6].x=oc.polar_center.x-L;
curraux[6].y=oc.polar_center.y+L;
curraux[6].z=oc.polar_center.z-L;

curraux_count=7;
}

else if(I2==1 && I3==1 && I5==1)//如果粒子与面2，面3和面5相交类似于角B1，补偿7个附加球
{

curraux[0].x=oc.polar_center.x-L;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y-L;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x-L;
curraux[2].y=oc.polar_center.y-L;
curraux[2].z=oc.polar_center.z;

curraux[3].x=oc.polar_center.x-L;
curraux[3].y=oc.polar_center.y-L;
curraux[3].z=oc.polar_center.z-L;

curraux[4].x=oc.polar_center.x;
curraux[4].y=oc.polar_center.y;
curraux[4].z=oc.polar_center.z-L;

curraux[5].x=oc.polar_center.x;
curraux[5].y=oc.polar_center.y-L;
curraux[5].z=oc.polar_center.z-L;

curraux[6].x=oc.polar_center.x-L;
curraux[6].y=oc.polar_center.y;
curraux[6].z=oc.polar_center.z-L;

curraux_count=7;
}

else if(I2==1 && I4==1 && I5==1)//如果粒子与面2，面4和面5相交类似于角C1，补偿7个附加球
{

curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z-L;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y-L;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x+L;
curraux[2].y=oc.polar_center.y;
curraux[2].z=oc.polar_center.z;

curraux[3].x=oc.polar_center.x+L;
curraux[3].y=oc.polar_center.y-L;
curraux[3].z=oc.polar_center.z;

curraux[4].x=oc.polar_center.x+L;
curraux[4].y=oc.polar_center.y;
curraux[4].z=oc.polar_center.z-L;

curraux[5].x=oc.polar_center.x;
curraux[5].y=oc.polar_center.y-L;
curraux[5].z=oc.polar_center.z-L;

curraux[6].x=oc.polar_center.x+L;
curraux[6].y=oc.polar_center.y-L;
curraux[6].z=oc.polar_center.z-L;

curraux_count=7;
}
else if(I1==1 && I4==1 && I5==1)//如果粒子与面1，面4和面5相交类似于角D1，补偿7个附加球
{

curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z-L;

curraux[1].x=oc.polar_center.x+L;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y+L;
curraux[2].z=oc.polar_center.z;

curraux[3].x=oc.polar_center.x+L;
curraux[3].y=oc.polar_center.y+L;
curraux[3].z=oc.polar_center.z;

curraux[4].x=oc.polar_center.x+L;
curraux[4].y=oc.polar_center.y;
curraux[4].z=oc.polar_center.z-L;

curraux[5].x=oc.polar_center.x;
curraux[5].y=oc.polar_center.y+L;
curraux[5].z=oc.polar_center.z-L;

curraux[6].x=oc.polar_center.x+L;
curraux[6].y=oc.polar_center.y+L;
curraux[6].z=oc.polar_center.z-L;

curraux_count=7;
}

else if(I1==1 && I3==1 && I6==1)//如果粒子与面1，面3和面6相交类似于角E1，补偿7个附加球
{

curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z+L;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y+L;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x-L;
curraux[2].y=oc.polar_center.y;
curraux[2].z=oc.polar_center.z;

curraux[3].x=oc.polar_center.x-L;
curraux[3].y=oc.polar_center.y+L;
curraux[3].z=oc.polar_center.z;

curraux[4].x=oc.polar_center.x;
curraux[4].y=oc.polar_center.y+L;
curraux[4].z=oc.polar_center.z+L;

curraux[5].x=oc.polar_center.x-L;
curraux[5].y=oc.polar_center.y;
curraux[5].z=oc.polar_center.z+L;

curraux[6].x=oc.polar_center.x-L;
curraux[6].y=oc.polar_center.y+L;
curraux[6].z=oc.polar_center.z+L;

curraux_count=7;
}
else if(I2==1 && I3==1 && I6==1)//如果粒子与面2，面3和面6相交类似于角F1，补偿7个附加球
{

curraux[0].x=oc.polar_center.x-L;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z;

curraux[1].x=oc.polar_center.x;
curraux[1].y=oc.polar_center.y-L;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y;
curraux[2].z=oc.polar_center.z+L;

curraux[3].x=oc.polar_center.x-L;
curraux[3].y=oc.polar_center.y-L;
curraux[3].z=oc.polar_center.z;

curraux[4].x=oc.polar_center.x-L;
curraux[4].y=oc.polar_center.y;
curraux[4].z=oc.polar_center.z+L;

curraux[5].x=oc.polar_center.x;
curraux[5].y=oc.polar_center.y-L;
curraux[5].z=oc.polar_center.z+L;

curraux[6].x=oc.polar_center.x-L;
curraux[6].y=oc.polar_center.y-L;
curraux[6].z=oc.polar_center.z+L;
curraux_count=7;
}
else if(I2==1 && I4==1 && I6==1)//如果粒子与面2，面4和面6相交类似于角G1，补偿7个附加球
{

curraux[0].x=oc.polar_center.x;
curraux[0].y=oc.polar_center.y;
curraux[0].z=oc.polar_center.z+L;

curraux[1].x=oc.polar_center.x+L;
curraux[1].y=oc.polar_center.y;
curraux[1].z=oc.polar_center.z;

curraux[2].x=oc.polar_center.x;
curraux[2].y=oc.polar_center.y-L;
curraux[2].z=oc.polar_center.z;

curraux[3].x=oc.polar_center.x;
curraux[3].y=oc.polar_center.y-L;
curraux[3].z=oc.polar_center.z+L;

curraux[4].x=oc.polar_center.x+L;
curraux[4].y=oc.polar_center.y;
curraux[4].z=oc.polar_center.z+L;

curraux[5].x=oc.polar_center.x+L;
curraux[5].y=oc.polar_center.y-L;
curraux[5].z=oc.polar_center.z;

curraux[6].x=oc.polar_center.x+L;
curraux[6].y=oc.polar_center.y-L;
curraux[6].z=oc.polar_center.z+L;

curraux_count=7;
}
for(int k=0;k<curraux_count;k++)
{
AuxB[aux_count].polar_center.x=curraux[k].x;
AuxB[aux_count].polar_center.y=curraux[k].y;
AuxB[aux_count].polar_center.z=curraux[k].z;
AuxB[aux_count].a_fa=oc.a_fa;
AuxB[aux_count].bei_ta=oc.bei_ta;
AuxB[aux_count].ga_ma=oc.ga_ma;
for(int i=0;i<N;i++)
{
AuxB[aux_count].polar_angle[i]=oc.polar_angle[i];
AuxB[aux_count].polar_radius[i]=oc.polar_radius[i];
}
GetTheLocatedCoordinate(AuxB[aux_count]);
extrapolate_the_polygon_to_3_D(AuxB[aux_count]);
DividePolygonToTriangles(AuxB[aux_count]);
aux_count++;
}
delete[] curraux;
}
bool PolygonsPercalation::PercolationDectect(double a01, double b01, double c01, double d01, polygon* C, polygon* D, int& Ccount, int& Dcount, int& AUXBcount, bool& correspondvalplane)
{
bool validplaneparticle=false; //bool型变量，申明当前考察面的对应的面无效，即validB中没有粒子与之相交
while (1)
{
bool particlesintersect=false; //bool型标示变量，表示在D中已经找不到粒子与C中粒子相交
for (int i=0; i<Dcount; i++)
{
if (D[i].polar_center.x==DELFLAG && D[i].polar_center.y==DELFLAG
&& D[i].polar_center.z==DELFLAG)//表示D[i]已经从D中删除,那么就不需要考察D[i],进行到下一个D[i+1]
continue;
for (int j=0; j<Ccount; j++)
{
if (IntersectParticlesJudegmentInaccurately(D[i], C[j])==1
&&IntersectParticlesJudegmentInaccurately(C[j], D[i])==1)
{
if(IntersectParticlesJudegementAccurately(D[i],C[j])==true)
{
particlesintersect=true; //表明在D中存在粒子D[i]与C中粒子相交
      C[Ccount]=D[i]; //并且将D[i]成员添加到C中
      Ccount++;
      D[i].polar_center.x=DELFLAG; //将D[i]做上标记，表示从D中将此成员删除
      D[i].polar_center.y=DELFLAG;
      D[i].polar_center.z=DELFLAG;
      break;
}
}
}
}//for (int i=0; i<Dcount; i++)
if (particlesintersect==false)
break;
}//while (1)

for (int m=0; m<Ccount; m++)//检测validB类中是否有粒子与当前面的对应面相交，若有,则达到percolation
{
if (IntersectParticleToPlane(C[m], a01, b01, c01, d01)==true)
{
correspondvalplane=true; //表示在validB类中有粒子与当前面的对应面相交
for (int k=0; k<AUXBcount; k++)//判断C类与考察面的对应面相交的粒子是否与附加粒子相交
{
if(IntersectParticlesJudegmentInaccurately(C[m], AuxB[k])==1
&&IntersectParticlesJudegmentInaccurately(AuxB[k],C[m])==1)
{
if(IntersectParticlesJudegementAccurately(C[m],AuxB[k]))
{
   validplaneparticle=true;
   break;
}
}
}
if(validplaneparticle==true)
break;
}
}
return validplaneparticle;
}
bool PolygonsPercalation::IntersectParticlesJudegmentInaccurately(polygon& P1,polygon& P2)//判断两个粒子是否相交的粗糙判定
{
double distance[3];
////////////////////////////////////////////////////////////////////////////////判定两粒子是否相交
for(int j=0;j<N;j++)
{
for(int k=0;k<3;k++)
{
distance[k]=P2.DivideTriangles[j].vertexs[k].x*P1.NoramalVector.x+P2.DivideTriangles[j].vertexs[k].y*P1.NoramalVector.y+P2.DivideTriangles[j].vertexs[k].z*P1.NoramalVector.z
-P1.polar_center.x*P1.NoramalVector.x-P1.polar_center.y*P1.NoramalVector.y-P1.polar_center.z*P1.NoramalVector.z;
}
if(distance[0]>0&&distance[1]>0&&distance[2]>0)
continue;
      if(distance[0]<0&&distance[1]<0&&distance[2]<0)
continue;
else //相交和共面情况返回为1
return 1;
}
return 0;//多边形中没有三角形相交，返回0
}
bool PolygonsPercalation::IntersectParticlesJudegementAccurately(polygon& P1,polygon& P2)
{
point IntersectingLineDirectionVector;//两个平面交线的方向向量
point PointInline;//交线上某点坐标
point PointOfInstersection[4];//两平面交线与两个三角形边的交点
double a,b,d1,d2;//参数a和b用来计算交线上的某点坐标,d1和d2是两三角形所在平面方程的常数项
double min_x1[2],max_x1[2];//用来确定大小的量
double t[4];//三角形与交线的标量值
double ModulusOfDirectionVector;
double distance[3];//顶点到平面的有向距离
double p[3];
double c;//变量c用来存储n1xn2，计算两平面法向量点成
double max_tempt_x=0,min_tempt_x=0,max_tempt_y=0,min_tempt_y=0;//用于比较大小的参数
double max_x,min_x,min_y,max_y;//在坐标面XOY投影的x和y的投影的最小和最大值
double MoudulusN1,MoudulusN2;//计算平面法向量N1和平面法向量N2的模的平方
////////////////////////////////////////////////////////////////////////////////计算两平面交线的方向向量
IntersectingLineDirectionVector.x=P1.NoramalVector.y*P2.NoramalVector.z-P1.NoramalVector.z*P2.NoramalVector.y;
IntersectingLineDirectionVector.y=P1.NoramalVector.z*P2.NoramalVector.x-P1.NoramalVector.x*P2.NoramalVector.z;
IntersectingLineDirectionVector.z=P1.NoramalVector.x*P2.NoramalVector.y-P1.NoramalVector.y*P2.NoramalVector.x;
   ModulusOfDirectionVector=sqrt(pow(IntersectingLineDirectionVector.x,2)+pow(IntersectingLineDirectionVector.y,2)+pow(IntersectingLineDirectionVector.z,2));
IntersectingLineDirectionVector.x/=ModulusOfDirectionVector;
IntersectingLineDirectionVector.y/=ModulusOfDirectionVector;
IntersectingLineDirectionVector.z/=ModulusOfDirectionVector;
////////////////////////////////////////////////////////////////////////////////计算两平面交线上任意一点
d1=P1.polar_center.x*P1.NoramalVector.x+P1.polar_center.y*P1.NoramalVector.y+P1.polar_center.z*P1.NoramalVector.z;
d2=P2.polar_center.x*P2.NoramalVector.x+P2.polar_center.y*P2.NoramalVector.y+P2.polar_center.z*P2.NoramalVector.z;
c=P1.NoramalVector.x*P2.NoramalVector.x+P1.NoramalVector.y+P2.NoramalVector.y+P1.NoramalVector.z*P2.NoramalVector.z;//计算两平面法向量点乘
MoudulusN1=pow(P1.NoramalVector.x,2)+pow(P1.NoramalVector.y,2)+pow(P1.NoramalVector.z,2);//计算P1平面法向量的模的平方
MoudulusN2=pow(P2.NoramalVector.x,2)+pow(P2.NoramalVector.y,2)+pow(P2.NoramalVector.z,2);//计算平面N2平面法向量的模的平方
a=(d2*c-d1*MoudulusN2)/(pow(c,2)-MoudulusN1*MoudulusN2);
b=(d1*c-d2*MoudulusN1)/(pow(c,2)-MoudulusN1*MoudulusN2);
PointInline.x=a*P1.NoramalVector.x+b*P2.NoramalVector.x;//交线上一点的坐标
PointInline.y=a*P1.NoramalVector.y+b*P2.NoramalVector.y;
PointInline.z=a*P1.NoramalVector.z+b*P2.NoramalVector.z;
for(int i=0;i<N;i++)
{
for(int j=0;j<N;j++)
{
for(int k=0;k<3;k++)
{
distance[k]=P2.DivideTriangles[j].vertexs[k].x*P1.NoramalVector.x+P2.DivideTriangles[j].vertexs[k].y*P1.NoramalVector.y+P2.DivideTriangles[j].vertexs[k].z*P1.NoramalVector.z
-P1.polar_center.x*P1.NoramalVector.x-P1.polar_center.y*P1.NoramalVector.y-P1.polar_center.z*P1.NoramalVector.z;
p[k]=P2.DivideTriangles[j].vertexs[k].x*IntersectingLineDirectionVector.x+P2.DivideTriangles[j].vertexs[k].y*IntersectingLineDirectionVector.y+
      P2.DivideTriangles[j].vertexs[k].z*IntersectingLineDirectionVector.z;
      }
//********************************************************当两个三角形共面时
if(distance[0]==0&&distance[1]==0&&distance[2]==0)
{
for(int t=0;t<3;t++)
{
P1.DivideTriangles[i].ProjectionVertexs[t].x=P1.DivideTriangles[i].vertexs[t].x;
P1.DivideTriangles[i].ProjectionVertexs[t].y=P1.DivideTriangles[i].vertexs[t].y;
P1.DivideTriangles[i].ProjectionVertexs[t].z=0;

P2.DivideTriangles[j].ProjectionVertexs[t].x=P2.DivideTriangles[j].vertexs[t].x;
P2.DivideTriangles[j].ProjectionVertexs[t].y=P2.DivideTriangles[j].vertexs[t].y;
P2.DivideTriangles[j].ProjectionVertexs[t].z=0;
}
for(int t=0;t<3;t++)
         {
if(min_tempt_x<P2.DivideTriangles[j].ProjectionVertexs[t].x)
            {max_tempt_x=P2.DivideTriangles[j].ProjectionVertexs[t].x;}
         else if(min_tempt_x>P2.DivideTriangles[j].ProjectionVertexs[t].x)
            {min_tempt_x=P2.DivideTriangles[j].ProjectionVertexs[t].x;}
         if(min_tempt_y<P2.DivideTriangles[j].ProjectionVertexs[t].y)
            {max_tempt_y=P2.DivideTriangles[j].ProjectionVertexs[t].y;}
         else if(min_tempt_y>P2.DivideTriangles[j].ProjectionVertexs[t].y)
            {min_tempt_y=P2.DivideTriangles[j].ProjectionVertexs[t].y;}
         }
min_x=min_tempt_x;max_x=max_tempt_x;min_y=min_tempt_y;max_y=max_tempt_y;
         for(int t=0;t<3;t++)
{
if(min_x<P1.DivideTriangles[i].ProjectionVertexs[t].x&&max_x>P1.DivideTriangles[i].ProjectionVertexs[t].x&&min_y<P1.DivideTriangles[i].ProjectionVertexs[t].y&&max_y>P1.DivideTriangles[i].ProjectionVertexs[t].y)
            {return(true);}//两共面三角形相交，返回1
         if(min_x==P1.DivideTriangles[i].ProjectionVertexs[t].x&&max_x==P1.DivideTriangles[i].ProjectionVertexs[t].x&&min_y==P1.DivideTriangles[i].ProjectionVertexs[t].y&&max_y==P1.DivideTriangles[i].ProjectionVertexs[t].y)
            {return(true);}//两共面三角形相交，返回1
}
}
if(distance[0]>0&&distance[1]<0&&distance[2]<0)//此时两平面交线与三角形边01，02边相交
{
t[0]=p[1]+(p[0]-p[1])*distance[1]/(distance[1]-distance[0]);
t[1]=p[2]+(p[0]-p[2])*distance[2]/(distance[2]-distance[0]);
}
if(distance[0]<0&&distance[1]>0&&distance[2]>0)//此时两平面交线与三角形边01，02边相交
{
t[0]=p[1]+(p[0]-p[1])*distance[1]/(distance[1]-distance[0]);
t[1]=p[2]+(p[0]-p[2])*distance[2]/(distance[2]-distance[0]);
}
if(distance[0]<0&&distance[1]<0&&distance[2]>0)//此时两平面交线与三角形边02，12边相交
{
t[0]=p[0]+(p[2]-p[0])*distance[0]/(distance[0]-distance[2]);
t[1]=p[1]+(p[2]-p[1])*distance[1]/(distance[1]-distance[2]);
}
if(distance[0]>0&&distance[1]>0&&distance[2]<0)//此时两平面交线与三角形边02，12边相交
{
t[0]=p[0]+(p[2]-p[0])*distance[0]/(distance[0]-distance[2]);
t[1]=p[1]+(p[2]-p[1])*distance[1]/(distance[1]-distance[2]);
}
if(distance[0]>0&&distance[2]>0&&distance[1]<0)//此时两平面交线与三角形01，12边相交
{
t[0]=p[0]+(p[1]-p[0])*distance[0]/(distance[0]-distance[1]);
t[1]=p[2]+(p[1]-p[2])*distance[2]/(distance[2]-distance[1]);
}
if(distance[0]<0&&distance[2]<0&&distance[1]>0)//此时两平面交线与三角形01，12边相交
{
t[0]=p[0]+(p[1]-p[0])*distance[0]/(distance[0]-distance[1]);
t[1]=p[2]+(p[1]-p[2])*distance[2]/(distance[2]-distance[1]);
}
////////////////////////////////////////////////////////////////////////////////////判定多边形P1与两平面的标量值
for(int k=0;k<3;k++)
{
distance[k]=P1.DivideTriangles[i].vertexs[k].x*P2.NoramalVector.x+P1.DivideTriangles[i].vertexs[k].y*P2.NoramalVector.y+P1.DivideTriangles[i].vertexs[k].z*P2.NoramalVector.z
-P2.polar_center.x*P2.NoramalVector.x-P2.polar_center.y*P2.NoramalVector.y-P2.polar_center.z*P2.NoramalVector.z;
p[k]=P1.DivideTriangles[i].vertexs[k].x*IntersectingLineDirectionVector.x+P1.DivideTriangles[i].vertexs[k].y*IntersectingLineDirectionVector.y+
      P1.DivideTriangles[i].vertexs[k].z*IntersectingLineDirectionVector.z;
      }
if(distance[0]>0&&distance[1]<0&&distance[2]<0)//此时两平面交线与三角形边01，02边相交
{
t[2]=p[1]+(p[0]-p[1])*distance[1]/(distance[1]-distance[0]);
t[3]=p[2]+(p[0]-p[2])*distance[2]/(distance[2]-distance[0]);
}
if(distance[0]<0&&distance[1]>0&&distance[2]>0)//此时两平面交线与三角形边01，02边相交
{
t[2]=p[1]+(p[0]-p[1])*distance[1]/(distance[1]-distance[0]);
t[3]=p[2]+(p[0]-p[2])*distance[2]/(distance[2]-distance[0]);
}
if(distance[0]<0&&distance[1]<0&&distance[2]>0)//此时两平面交线与三角形边02，12边相交
{
t[2]=p[0]+(p[2]-p[0])*distance[0]/(distance[0]-distance[2]);
t[3]=p[1]+(p[2]-p[1])*distance[1]/(distance[1]-distance[2]);
}
if(distance[0]>0&&distance[1]>0&&distance[2]<0)//此时两平面交线与三角形边02，12边相交
{
t[2]=p[0]+(p[2]-p[0])*distance[0]/(distance[0]-distance[2]);
t[3]=p[1]+(p[2]-p[1])*distance[1]/(distance[1]-distance[2]);
}
if(distance[0]>0&&distance[2]>0&&distance[1]<0)//此时两平面交线与三角形01，12边相交
{
t[2]=p[0]+(p[1]-p[0])*distance[0]/(distance[0]-distance[1]);
t[3]=p[2]+(p[1]-p[2])*distance[2]/(distance[2]-distance[1]);
}
if(distance[0]<0&&distance[2]<0&&distance[1]>0)//此时两平面交线与三角形01，12边相交
{
t[2]=p[0]+(p[1]-p[0])*distance[0]/(distance[0]-distance[1]);
t[3]=p[2]+(p[1]-p[2])*distance[2]/(distance[2]-distance[1]);
}
for(int r=0;r<4;r++)//计算两平面的交线与两三角形边的交点
{
PointOfInstersection[r].x=t[r]*IntersectingLineDirectionVector.x+PointInline.x;
PointOfInstersection[r].y=t[r]*IntersectingLineDirectionVector.y+PointInline.y;
PointOfInstersection[r].z=t[r]*IntersectingLineDirectionVector.z+PointInline.z;
}
if(PointOfInstersection[0].x>PointOfInstersection[1].x)
{
min_x1[0]=PointOfInstersection[1].x;
max_x1[0]=PointOfInstersection[0].x;
}
if(PointOfInstersection[0].x<PointOfInstersection[1].x)
{
min_x1[0]=PointOfInstersection[0].x;
max_x1[0]=PointOfInstersection[1].x;
}
if(PointOfInstersection[2].x>PointOfInstersection[3].x)
{
min_x1[1]=PointOfInstersection[3].x;
max_x1[1]=PointOfInstersection[2].x;
}
if(PointOfInstersection[2].x<PointOfInstersection[3].x)
{
min_x1[1]=PointOfInstersection[2].x;
max_x1[1]=PointOfInstersection[3].x;
}
if(min_x1[0]<=max_x1[1]&&max_x1[0]>=max_x1[1])
{
return true;//此时两不共面三角形相交，返回true
}
else if(min_x1[0]<=min_x1[1]&&min_x1[1]<=max_x1[0])
{
return true;//此时两不共面三角形相交，返回true
}
else
continue;
}
}
return false;
}
int main()
{
srand((int)time(NULL));
clock_t start, end; //程序执行时间(ms)
double TotalTime;
start=clock();
PolygonsPercalation P;
P.GetDate();
P.PercolationSimulation();
end=clock();
TotalTime=((end-start)*0.001)/3600;
cout<<"TotalTime="<<TotalTime<<" hours"<<endl;
return 0;
}

账号		自动登录	找回密码
密码			立即注册

[技术交流] 出现无法读取内存情况怎么办？

马上注册，结交更多好友，享用更多功能^_^