C++ HCDT类代码示例

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_datac_set_ri(HADT d, HCDT re, HCDT im)
{
	long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
	d->Create(nx,ny,nz);
#pragma omp parallel for
	for(long i=0;i<nx*ny*nz;i++)	d->a[i] = dual(re->vthr(i),im->vthr(i));
}

HADT MGL_EXPORT mgl_datac_sum(HCDT dat, const char *dir)
{
	if(!dir || *dir==0)	return 0;
	long nx=dat->GetNx(),ny=dat->GetNy(),nz=dat->GetNz();
	long p[3]={nx,ny,nz};
	dual *b = new dual[nx*ny*nz];
	dual *c = new dual[nx*ny*nz];

	const mglDataC *d=dynamic_cast<const mglDataC *>(dat);
	if(d)	memcpy(c,d->a,nx*ny*nz*sizeof(dual));
	else
#pragma omp parallel for
		for(long i=0;i<nx*ny*nz;i++)	c[i]=dat->vthr(i);

	if(strchr(dir,'z') && nz>1)
	{
		mglStartThreadC(mgl_sumc_z,0,nx*ny,b,c,0,p);
		memcpy(c,b,nx*ny*sizeof(mreal));	p[2] = 1;
	}
	if(strchr(dir,'y') && ny>1)
	{
		mglStartThreadC(mgl_sumc_y,0,nx*p[2],b,c,0,p);
		memcpy(c,b,nx*p[2]*sizeof(mreal));	p[1] = p[2];	p[2] = 1;
	}
	if(strchr(dir,'x') && nx>1)
	{
		mglStartThreadC(mgl_sumc_x,0,p[1]*p[2],b,c,0,p);
		p[0] = p[1];	p[1] = p[2];	p[2] = 1;
	}
	mglDataC *r=new mglDataC(p[0],p[1],p[2]);
	memcpy(r->a,b,p[0]*p[1]*p[2]*sizeof(dual));
	delete []b;	delete []c;	return r;
}

//-----------------------------------------------------------------------------
HADT MGL_EXPORT mgl_datac_column(HCDT dat, const char *eq)
{
	const mglData *dd=dynamic_cast<const mglData *>(dat);
	std::vector<mglDataA*> list;
	if(dd && dd->id.length()>0)	for(size_t i=0;i<dd->id.length();i++)
	{
		mglDataT *col = new mglDataT(*dat);
		col->SetInd(i,dd->id[i]);
		list.push_back(col);
	}
	const mglDataC *dc=dynamic_cast<const mglDataC *>(dat);
	if(dc && dc->id.length()>0)	for(size_t i=0;i<dc->id.length();i++)
	{
		mglDataT *col = new mglDataT(*dat);
		col->SetInd(i,dc->id[i]);
		list.push_back(col);
	}
	if(list.size()==0)	return 0;	// no named columns
	mglDataV *t = new mglDataV(dat->GetNy(),dat->GetNz());
	t->s=L"#$mgl";	list.push_back(t);
	mglDataC *r = new mglDataC;
	r->Set(mglFormulaCalcC(eq,list));
	for(size_t i=0;i<list.size();i++)	delete list[i];
	return r;
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_data_grid_xy(HMDT d, HCDT xdat, HCDT ydat, HCDT zdat, mreal x1, mreal x2, mreal y1, mreal y2)
{ // NOTE: only for mglData
	const mglData *x = dynamic_cast<const mglData *>(xdat);
	const mglData *y = dynamic_cast<const mglData *>(ydat);
	long n=xdat->GetNN();
	if((n<3) || (ydat->GetNN()!=n) || (zdat->GetNN()!=n))	return;

	mglData *nums = mgl_triangulation_2d(xdat,ydat);
	if(nums->nx<3)	{	delete nums;	return;	}
	long nn = nums->ny, par[3]={d->nx,d->ny,d->nz};
	mreal xx[4]={x1,(d->nx-1.)/(x2-x1), y1,(d->ny-1.)/(y2-y1)};

	mreal *xc=new mreal[n], *yc=new mreal[n];
	if(x && y)
#pragma omp parallel for
		for(long i=0;i<n;i++)
		{	xc[i]=xx[1]*(x->a[i]-xx[0]);	yc[i]=xx[3]*(y->a[i]-xx[2]);	}
	else
#pragma omp parallel for
		for(long i=0;i<n;i++)
		{	xc[i]=xx[1]*(xdat->vthr(i)-xx[0]);	yc[i]=xx[3]*(ydat->vthr(i)-xx[2]);	}
	long tmp = d->nx*d->ny*d->nz;
#pragma omp parallel for
	for(long i=0;i<tmp;i++) d->a[i] = NAN;

	mglStartThread(mgl_grid_t,0,nn,d->a,xc,yc,par,zdat,nums->a);
	delete nums;	delete []xc;	delete []yc;
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_cont_z_val(HMGL gr, HCDT v, HCDT a, const char *sch, double sv, const char *opt)
{
	long n=a->GetNx(),m=a->GetNy();
	if(n<2 || m<2)	{	gr->SetWarn(mglWarnLow,"ContZ");	return;	}
	gr->SaveState(opt);
	if(mgl_isnan(sv))	sv = gr->GetOrgZ('z');
	if(sv<gr->Min.z || sv>gr->Max.z)	{	gr->SetWarn(mglWarnSlc,"ContZ");	gr->LoadState();	return;	}
	static int cgid=1;	gr->StartGroup("ContZ",cgid++);
	mglDataV xx,yy,zz;	mglData aa;

	int text=0;
	if(mglchr(sch,'t'))	text=1;
	if(mglchr(sch,'T'))	text=2;
	long ss=gr->AddTexture(sch);
	gr->SetPenPal(sch);

	a = fill_slice_z(gr,sv,a,xx,yy,zz,aa);
#pragma omp parallel for
	for(long i=0;i<v->GetNx();i++)
	{
		mreal v0 = v->v(i);
		mgl_cont_gen(gr,v0,a,&xx,&yy,&zz,gr->GetC(ss,v0),text,0);
	}
	gr->EndGroup();
}

//-----------------------------------------------------------------------------
//
//	Crust series
//
//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_crust(HMGL gr, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
{
	if(y->GetNx()!=x->GetNx() || z->GetNx()!=x->GetNx())
	{	gr->SetWarn(mglWarnDim,"Crust");	return;	}
	HMDT nums = mgl_triangulation_3d(x, y, z);
	mgl_triplot_xyzc(gr,nums,x,y,z,z,sch,opt);
	mgl_delete_data(nums);
}

//-----------------------------------------------------------------------------
HMDT MGL_EXPORT mgl_datac_arg(HCDT d)
{
	long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
	mglData *r=new mglData(nx,ny,nz);
	const mglDataC *dd = dynamic_cast<const mglDataC*>(d);
	if(dd)
#pragma omp parallel for
		for(long i=0;i<nx*ny*nz;i++)	r->a[i] = arg(dd->a[i]);
	return r;
}

//-----------------------------------------------------------------------------
//
//	TriContV series
//
//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_tricontv_xyzcv(HMGL gr, HCDT v, HCDT nums, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, const char *opt)
{
	mglDataV zz(x->GetNx(),x->GetNy());
	if(!z)	z = &zz;
	if(mgl_check_trig(gr,nums,x,y,z,a,"TriContV"))	return;

	gr->SaveState(opt);
	static int cgid=1;	gr->StartGroup("TriContV",cgid++);
	bool fixed=(mglchr(sch,'_')) || (gr->Min.z==gr->Max.z);
	long s=gr->AddTexture(sch);
	gr->SetPenPal(sch);

	for(long k=0;k<v->GetNx();k++)
	{
		mreal v0 = v->v(k);		zz.Fill(fixed ? gr->Min.z : v0);
		mreal dv = (gr->Max.c-gr->Min.c)/8, c = gr->GetC(s,v0);
		if(k>0)	dv = v->v(k-1)-v->v(k);
		else if(k<v->GetNx()-1)	dv = v->v(k)-v->v(k+1);
		if(fixed)	dv=-dv;

		const std::vector<mglSegment> curvs = mgl_get_curvs(gr,mgl_tri_lines(v0,nums,a,x,y,fixed?&zz:z));
		for(size_t i=0;i<curvs.size();i++)
		{
			const std::list<mglPoint> &pp=curvs[i].pp;
			long f2=-1,g2=-1;
			for(std::list<mglPoint>::const_iterator it=pp.begin(); it != pp.end(); ++it)
			{
				mglPoint p=*it,q(p.y,-p.x);
				long f1 = f2;	f2 = gr->AddPnt(p,c,q);	p.z+=dv;
				long g1 = g2;	g2 = gr->AddPnt(p,c,q);
				gr->quad_plot(f1,g1,f2,g2);
			}
		}
	}
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_datac_set_ap(HADT d, HCDT a, HCDT p)
{
	long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
	d->Create(nx,ny,nz);
#pragma omp parallel for
	for(long i=0;i<nx*ny*nz;i++)
	{
		register mreal aa=a->vthr(i), pp=p->vthr(i);
		d->a[i] = dual(aa*cos(pp), aa*sin(pp));
	}
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_dens_z(HMGL gr, HCDT a, const char *sch, double sv, const char *opt)
{
	long n=a->GetNx(),m=a->GetNy();
	if(n<2 || m<2)	{	gr->SetWarn(mglWarnLow,"DensZ");	return;	}
	gr->SaveState(opt);
	if(mgl_isnan(sv))	sv = gr->GetOrgZ('z');
	if(sv<gr->Min.z || sv>gr->Max.z)	{	gr->SetWarn(mglWarnSlc,"DensZ");	gr->LoadState();	return;	}
	mglDataV xx,yy,zz;	mglData aa;
	a = fill_slice_z(gr,sv,a,xx,yy,zz,aa);
	mgl_surf_gen(gr, &xx,&yy,&zz,a, 0, sch);
}

//-----------------------------------------------------------------------------
//
//	IFS series
//
//-----------------------------------------------------------------------------
void MGL_NO_EXPORT mgl_ifs_2d_point(HCDT A, mreal& x, mreal& y, mreal amax)
{
	long i, n=A->GetNy();
	mreal r = amax*mgl_rnd(), sum_prob = 0, x1;
	for(i=0;i<n;i++)
	{
		sum_prob += A->v(6,i);
		if(r<sum_prob)	break;
	}
	x1= A->v(0,i)*x + A->v(1,i)*y + A->v(4,i);
	y = A->v(2,i)*x + A->v(3,i)*y + A->v(5,i);	x = x1;
}

//-----------------------------------------------------------------------------
void MGL_NO_EXPORT mgl_ifs_3d_point(HCDT A, mreal& x, mreal& y, mreal& z, mreal amax)
{
	int i, n=A->GetNy();
	mreal r = amax*mgl_rnd(), sum_prob = 0, x1, y1;
	for (i=0; i<n; i++)
	{
		sum_prob += A->v(12,i);
		if(r < sum_prob)  break;
	}
	x1= A->v(0,i)*x + A->v(1,i)*y + A->v(2,i)*z + A->v(9,i);
	y1= A->v(3,i)*x + A->v(4,i)*y + A->v(5,i)*z + A->v(10,i);
	z = A->v(6,i)*x + A->v(7,i)*y + A->v(8,i)*z + A->v(11,i);
	x = x1;	y = y1;
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_datac_set(HADT d, HCDT a)
{
	if(!a)	return;
	const mglDataC *dd = dynamic_cast<const mglDataC *>(a);	// faster for mglData
	mgl_datac_create(d, a->GetNx(), a->GetNy(), a->GetNz());
	if(dd)	// this one should be much faster
		memcpy(d->a, dd->a, d->nx*d->ny*d->nz*sizeof(dual));
	else	// very inefficient!!!
	{
		register long i,j,k;
		for(k=0;k<d->nz;k++)	for(j=0;j<d->ny;j++)	for(i=0;i<d->nx;i++)
			d->a[i+d->nx*(j+d->ny*k)] = a->v(i,j,k);
	}
}

void MGL_EXPORT mgl_datac_modify_vw(HADT d, const char *eq,HCDT vdat,HCDT wdat)
{
	const mglDataC *v = dynamic_cast<const mglDataC *>(vdat);
	const mglDataC *w = dynamic_cast<const mglDataC *>(wdat);
	long nn = d->nx*d->ny*d->nz, par[3]={d->nx,d->ny,d->nz};
	if(vdat && vdat->GetNN()!=nn)	return;
	if(wdat && wdat->GetNN()!=nn)	return;
	mglFormulaC f(eq);
	if(v && w)	mglStartThreadC(mgl_cmodify,0,nn,d->a,v->a,w->a,par,&f);
	else if(vdat && wdat)	mglStartThreadV(mgl_cmodify_gen,nn,d->a,vdat,wdat,par,&f);
	else if(v)	mglStartThreadC(mgl_cmodify,0,nn,d->a,v->a,0,par,&f);
	else if(vdat)	mglStartThreadV(mgl_cmodify_gen,nn,d->a,vdat,0,par,&f);
	else	mglStartThreadC(mgl_cmodify,0,nn,d->a,0,0,par,&f);
}

//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_triplot_xy(HMGL gr, HCDT nums, HCDT x, HCDT y, const char *sch, const char *opt)
{
	gr->SaveState(opt);
	mglData z(x->GetNx());
	mreal zm = gr->AdjustZMin();	z.Fill(zm,zm);
	mgl_triplot_xyzc(gr,nums,x,y,&z,&z,sch,0);
}