C++ HCDT::vthr方法代码示例

//-----------------------------------------------------------------------------
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;
}

//-----------------------------------------------------------------------------
HADT MGL_EXPORT mgl_datac_combine(HCDT d1, HCDT d2)
{
	long n1=d1->GetNy(),n2=d2->GetNx(),nx=d1->GetNx();
	if(d1->GetNz()>1 || (n1>1 && d2->GetNy()>1) || d2->GetNz()>1)	return 0;	// wrong dimensions
	mglDataC *r=new mglDataC;
	bool dim2=true;
	if(n1==1)	{	n1=n2;	n2=d2->GetNy();	dim2 = false;	}
	r->Create(nx,n1,n2);
	if(dim2)	n1*=nx;	else	{	n2*=n1;	n1=nx;	}

	const mglDataC *c1=dynamic_cast<const mglDataC *>(d1);
	const mglDataC *c2=dynamic_cast<const mglDataC *>(d2);
	if(c1 && c2)
#pragma omp parallel for collapse(2)
		for(long j=0;j<n2;j++)	for(long i=0;i<n1;i++)
			r->a[i+n1*j] = c1->a[i]*c2->a[j];
	else if(c1)
#pragma omp parallel for collapse(2)
		for(long j=0;j<n2;j++)	for(long i=0;i<n1;i++)
			r->a[i+n1*j] = c1->a[i]*d2->vthr(j);
	else if(c2)
#pragma omp parallel for collapse(2)
		for(long j=0;j<n2;j++)	for(long i=0;i<n1;i++)
			r->a[i+n1*j] = d1->vthr(i)*c2->a[j];
	else
#pragma omp parallel for collapse(2)
		for(long j=0;j<n2;j++)	for(long i=0;i<n1;i++)
			r->a[i+n1*j] = d1->vthr(i)*d2->vthr(j);
	return r;
}

//-----------------------------------------------------------------------------
//
//	Dots series
//
//-----------------------------------------------------------------------------
void MGL_EXPORT mgl_dots_ca(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT c, HCDT a, const char *sch, const char *opt)
{
	long n = x->GetNN(), d, k=1;
	if(x->GetNz()>1) 	k=3;		else if(x->GetNy()>1)	k=2;

	if(y->GetNN()!=n || z->GetNN()!=n || c->GetNN()!=n || (a && a->GetNN()!=n))
	{	gr->SetWarn(mglWarnDim,"Dots");	return;	}
	gr->SaveState(opt);

	d = gr->MeshNum>0 ? mgl_ipow(gr->MeshNum+1,k) : n;
	d = n>d ? n/d:1;

	static int cgid=1;	gr->StartGroup("Dots",cgid++);
	char mk=gr->SetPenPal(sch);
	long ss=gr->AddTexture(sch);
	if(mk==0)	mk='.';
	gr->Reserve(n);

	for(long i=0;i<n;i+=d)
	{
		mglPoint p = mglPoint(x->vthr(i),y->vthr(i),z->vthr(i));
		long pp = gr->AddPnt(p,gr->GetC(ss,c->vthr(i)),mglPoint(NAN),a?gr->GetA(a->vthr(i)):-1);
		gr->mark_plot(pp, mk);
	}
	gr->EndGroup();
}

//-----------------------------------------------------------------------------
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_evaluate(HCDT dat, HCDT idat, HCDT jdat, HCDT kdat, int norm)
{
	if(!idat || (jdat && jdat->GetNN()!=idat->GetNN()) || (kdat && kdat->GetNN()!=idat->GetNN()))	return 0;
	const mglData *dd=dynamic_cast<const mglData *>(dat);
	const mglDataC *dc=dynamic_cast<const mglDataC *>(dat);
	long nx=dat->GetNx(), ny=dat->GetNy(), nz=dat->GetNz();
	mglDataC *r=new mglDataC(idat->GetNx(),idat->GetNy(),idat->GetNz());
	if(dd)
#pragma omp parallel for
		for(long i=0;i<idat->GetNN();i++)
		{
			mreal x=idat->vthr(i), y=jdat?jdat->vthr(i):0, z=kdat?kdat->vthr(i):0;
			r->a[i] = mgl_isnum(x*y*z)?mglSpline3st<mreal>(dd->a,nx,ny,nz, x,y,z):NAN;
		}
	else if(dc)
#pragma omp parallel for
		for(long i=0;i<idat->GetNN();i++)
		{
			mreal x=idat->vthr(i), y=jdat?jdat->vthr(i):0, z=kdat?kdat->vthr(i):0;
			r->a[i] = mgl_isnum(x*y*z)?mglSpline3st<dual>(dc->a,nx,ny,nz, x,y,z):NAN;
		}
	else
#pragma omp parallel for
		for(long i=0;i<idat->GetNN();i++)
		{
			mreal x=idat->vthr(i), y=jdat?jdat->vthr(i):0, z=kdat?kdat->vthr(i):0;
			r->a[i] = mgl_isnum(x*y*z)?mgl_data_linear(dat, x,y,z):NAN;;
		}
	return r;
}

//-----------------------------------------------------------------------------
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));
	}
}

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;
}

//-----------------------------------------------------------------------------
//
//	DataGrid
//
//-----------------------------------------------------------------------------
MGL_NO_EXPORT void *mgl_grid_t(void *par)
{
	mglThreadD *t=(mglThreadD *)par;
	long nx=t->p[0],ny=t->p[1];
	mreal *b=t->a;
	const mreal *x=t->b, *y=t->c, *d=t->d;
	HCDT zdat = (HCDT) t->v;
#if !MGL_HAVE_PTHREAD
#pragma omp parallel for
#endif
	for(long i0=t->id;i0<t->n;i0+=mglNumThr)
	{
		register long k1 = long(d[3*i0]+0.5), k2 = long(d[3*i0+1]+0.5), k3 = long(d[3*i0+2]+0.5);
		mreal dxu,dxv,dyu,dyv;
		mreal z1=zdat->vthr(k1), z2=zdat->vthr(k2), z3=zdat->vthr(k3);
		mglPoint d1=mglPoint(x[k2]-x[k1],y[k2]-y[k1],z2-z1), d2=mglPoint(x[k3]-x[k1],y[k3]-y[k1],z3-z1), p;

		dxu = d2.x*d1.y - d1.x*d2.y;
		if(fabs(dxu)<1e-5) continue; // points lies on the same line
		dyv =-d1.x/dxu; dxv = d1.y/dxu;
		dyu = d2.x/dxu; dxu =-d2.y/dxu;

		long x1,y1,x2,y2;
		x1 = long(mgl_min(mgl_min(x[k1],x[k2]),x[k3])); // bounding box
		y1 = long(mgl_min(mgl_min(y[k1],y[k2]),y[k3]));
		x2 = long(mgl_max(mgl_max(x[k1],x[k2]),x[k3]));
		y2 = long(mgl_max(mgl_max(y[k1],y[k2]),y[k3]));
		x1 = x1>0 ? x1:0; x2 = x2<nx ? x2:nx-1;
		y1 = y1>0 ? y1:0; y2 = y2<ny ? y2:ny-1;
		if((x1>x2) | (y1>y2)) continue;

		register mreal u,v,xx,yy, x0 = x[k1], y0 = y[k1];
		register long i,j;
		for(i=x1;i<=x2;i++) for(j=y1;j<=y2;j++)
		{
			xx = (i-x0); yy = (j-y0);
			u = dxu*xx+dyu*yy; v = dxv*xx+dyv*yy;
			if((u<0) | (v<0) | (u+v>1)) continue;
			b[i+nx*j] = z1 + d1.z*u + d2.z*v;
		}
	}
	return 0;
}

HMDT MGL_EXPORT mgl_triangulation_2d(HCDT x, HCDT y)
{
	mglData *nums=0;
	long n = x->GetNN();
	if(y->GetNN()!=n)	return nums;
	// use s-hull here
	std::vector<Shx> pts;
	std::vector<long> out;
	Shx pt;

	double mx = 0, my = 0;
	for(long i=0;i<n;i++)
	{
		register double t;
		t = fabs(x->vthr(i));	if(t>mx)	mx=t;
		t = fabs(y->vthr(i));	if(t>my)	my=t;
	}
	mx *= 1e-15;	my *= 1e-15;
	for(long i=0;i<n;i++)
	{
		pt.r = x->vthr(i);	pt.c = y->vthr(i);
		if(mgl_isbad(pt.r) || mgl_isbad(pt.c))	continue;
		if(fabs(pt.r)<mx)	pt.r=0;
		if(fabs(pt.c)<my)	pt.c=0;
		pt.id = i;    pts.push_back(pt);
	}

	std::vector<Triad> triads;
	if(de_duplicate(pts, out))
		mgl_set_global_warn("There are duplicated points for triangulation.");
	s_hull_pro(pts, triads);
	long m = triads.size();
	nums=new mglData(3,m);
	for(long i=0;i<m;i++)
	{
		nums->a[3*i]   = triads[i].a;
		nums->a[3*i+1] = triads[i].b;
		nums->a[3*i+2] = triads[i].c;
	}
	return nums;
}

//-----------------------------------------------------------------------------
HMDT MGL_EXPORT mgl_datac_abs(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] = abs(dd->a[i]);
	else
#pragma omp parallel for
		for(long i=0;i<nx*ny*nz;i++)	r->a[i] = fabs(d->vthr(i));
	return r;
}

//-----------------------------------------------------------------------------
HADT MGL_EXPORT mgl_datac_subdata_ext(HCDT d, HCDT xx, HCDT yy, HCDT zz)
{
	if(!xx || !yy || !zz)
	{
		mglData tmp;	tmp.a[0]=-1;
		return mgl_datac_subdata_ext(d,xx?xx:&tmp,yy?yy:&tmp,zz?zz:&tmp);
	}
	
	long n=0,m=0,l=0,j,k;
	bool ix=false, iy=false, iz=false;
	if(xx->GetNz()>1)	// 3d data
	{
		n = xx->GetNx();	m = xx->GetNy();	l = xx->GetNz();
		j = yy->GetNN();	if(j>1 && j!=n*m*l)	return 0;	// wrong sizes
		k = zz->GetNN();	if(k>1 && k!=n*m*l)	return 0;	// wrong sizes
		ix = true;	iy = j>1;	iz = k>1;
	}
	else if(yy->GetNz()>1)
	{
		n = yy->GetNx();	m = yy->GetNy();	l = yy->GetNz();
		j = xx->GetNN();	if(j>1 && j!=n*m*l)	return 0;	// wrong sizes
		k = zz->GetNN();	if(k>1 && k!=n*m*l)	return 0;	// wrong sizes
		iy = true;	ix = j>1;	iz = k>1;
	}
	else if(zz->GetNz()>1)
	{
		n = zz->GetNx();	m = zz->GetNy();	l = zz->GetNz();
		j = yy->GetNN();	if(j>1 && j!=n*m*l)	return 0;	// wrong sizes
		k = xx->GetNN();	if(k>1 && k!=n*m*l)	return 0;	// wrong sizes
		iz = true;	iy = j>1;	ix = k>1;
	}
	else if(xx->GetNy()>1)	// 2d data
	{
		n = xx->GetNx();	m = xx->GetNy();	l = 1;
		j = yy->GetNx()*yy->GetNy();	if(j>1 && j!=n*m)	return 0;	// wrong sizes
		k = zz->GetNx()*zz->GetNy();	if(k>1 && k!=n*m)	return 0;	// wrong sizes
		ix = true;	iy = j>1;	iz = k>1;
	}
	else if(yy->GetNy()>1)
	{
		n = yy->GetNx();	m = yy->GetNy();	l = 1;
		j = xx->GetNx()*xx->GetNy();	if(j>1 && j!=n*m)	return 0;	// wrong sizes
		k = zz->GetNx()*zz->GetNy();	if(k>1 && k!=n*m)	return 0;	// wrong sizes
		iy = true;	ix = j>1;	iz = k>1;
	}
	else if(zz->GetNy()>1)
	{
		n = zz->GetNx();	m = zz->GetNy();	l = 1;
		j = yy->GetNx()*yy->GetNy();	if(j>1 && j!=n*m)	return 0;	// wrong sizes
		k = xx->GetNx()*xx->GetNy();	if(k>1 && k!=n*m)	return 0;	// wrong sizes
		iz = true;	iy = j>1;	ix = k>1;
	}
	long nx=d->GetNx(),ny=d->GetNy(),nz=d->GetNz();
	long vx=long(xx->v(0)), vy=long(yy->v(0)), vz=long(zz->v(0));
	const mglDataC *dd = dynamic_cast<const mglDataC *>(d);
	if(n*m*l>1)	// this is 2d or 3d data
	{
		mglDataV tx(n,m,l),ty(n,m,l),tz(n,m,l);
		if(!ix)	{	xx = &tx;	if(vx>=0)	tx.Fill(vx);	else tx.All();	}
		if(!iy)	{	yy = &ty;	if(vy>=0)	ty.Fill(vy);	else ty.All();	}
		if(!iz)	{	zz = &tz;	if(vz>=0)	tz.Fill(vz);	else tz.All();	}
		mglDataC *r=new mglDataC(n,m,l);
		if(dd)
#pragma omp parallel for
			for(long i0=0;i0<n*m*l;i0++)
			{
				register long x=long(0.5+xx->vthr(i0)), y=long(0.5+yy->vthr(i0)), z=long(0.5+zz->vthr(i0));
				r->a[i0] = (x>=0 && x<nx && y>=0 && y<ny && z>=0 && z<nz)?dd->a[x+nx*(y+ny*z)]:NAN;
			}
		else
#pragma omp parallel for
			for(long i0=0;i0<n*m*l;i0++)
			{
				register long x=long(0.5+xx->vthr(i0)), y=long(0.5+yy->vthr(i0)), z=long(0.5+zz->vthr(i0));
				r->a[i0] = (x>=0 && x<nx && y>=0 && y<ny && z>=0 && z<nz)?d->v(x,y,z):NAN;
			}
		return r;
	}
	// this is 1d data -> try as normal SubData()
	mglDataV tx(nx),ty(ny),tz(nz);	tx.Fill(0,nx-1);	ty.Fill(0,ny-1);	tz.Fill(0,nz-1);
	if(xx->GetNx()>1 || vx>=0)	n=xx->GetNx();	else	{	n=nx;	xx = &tx;	}
	if(yy->GetNx()>1 || vy>=0)	m=yy->GetNx();	else	{	m=ny;	yy = &ty;	}
	if(zz->GetNx()>1 || vz>=0)	l=zz->GetNx();	else	{	l=nz;	zz = &tz;	}
	mglDataC *r=new mglDataC(n,m,l);
	if(dd)
#pragma omp parallel for collapse(3)
		for(long k=0;k<l;k++)	for(long j=0;j<m;j++)	for(long i=0;i<n;i++)
		{
			register long x=long(0.5+xx->v(i)), y=long(0.5+yy->v(j)), z=long(0.5+zz->v(k));
			r->a[i+n*(j+m*k)] = (x>=0 && x<nx && y>=0 && y<ny && z>=0 && z<nz)?dd->a[x+nx*(y+ny*z)]:NAN;
		}
	else
#pragma omp parallel for collapse(3)
		for(long k=0;k<l;k++)	for(long j=0;j<m;j++)	for(long i=0;i<n;i++)
		{
			register long x=long(0.5+xx->v(i)), y=long(0.5+yy->v(j)), z=long(0.5+zz->v(k));
			r->a[i+n*(j+m*k)] = (x>=0 && x<nx && y>=0 && y<ny && z>=0 && z<nz)?d->v(x,y,z):NAN;
		}
	if(m==1)	{	r->ny=r->nz;	r->nz=1;	}// "squeeze" dimensions
//.........这里部分代码省略.........

//-----------------------------------------------------------------------------
HADT MGL_EXPORT mgl_datac_momentum(HCDT dat, char dir, const char *how)
{
	if(!how || !(*how) || !strchr("xyz",dir))	return 0;
	long nx=dat->GetNx(),ny=dat->GetNy(),nz=dat->GetNz();
	mglDataV x(nx,ny,nz, 0,1,'x');	x.s=L"x";
	mglDataV y(nx,ny,nz, 0,1,'y');	y.s=L"y";
	mglDataV z(nx,ny,nz, 0,1,'z');	z.s=L"z";
	mglDataC u(dat);	u.s=L"u";	// NOTE slow !!!
	std::vector<mglDataA*> list;
	list.push_back(&x);	list.push_back(&y);	list.push_back(&z);	list.push_back(&u);
	mglDataC res=mglFormulaCalcC(how,list);

	mglDataC *b=0;
	if(dir=='x')
	{
		b=new mglDataC(nx);
#pragma omp parallel for
		for(long i=0;i<nx;i++)
		{
			register dual i1=0,i0=0;
			for(long j=0;j<ny*nz;j++)
			{
				register dual u=dat->vthr(i+nx*j);
				i0 += u;	i1 += u*res.a[i+nx*j];
			}
			b->a[i] = i0!=mreal(0) ? i1/i0 : 0;
		}
	}
	if(dir=='y')
	{
		b=new mglDataC(ny);
#pragma omp parallel for
		for(long i=0;i<ny;i++)
		{
			register dual i1=0,i0=0;
			for(long k=0;k<nz;k++)	for(long j=0;j<nx;j++)
			{
				register dual u=dat->v(j,i,k);
				i0 += u;	i1 += u*res.a[j+nx*(i+ny*k)];
			}
			b->a[i] = i0!=mreal(0) ? i1/i0 : 0;
		}
	}
	if(dir=='z')
	{
		long nn=nx*ny;
		b=new mglDataC(nz);
#pragma omp parallel for
		for(long i=0;i<nz;i++)
		{
			register dual i1=0,i0=0;
			for(long j=0;j<nn;j++)
			{
				register dual u=dat->vthr(j+nn*i);
				i0 += u;	i1 += u*res.a[j+nn*i];
			}
			b->a[i] = i0!=mreal(0) ? i1/i0 : 0;
		}
	}
	return b;
}