24小时热门版块排行榜

返回列表

焊接网格划分

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 906.1
帖子: 24
在线: 39.2小时
虫号: 2111966

[交流] 激光反射吸收热源模型！！！！代码重量爆料！！！！！求解答！！！！！！！！！

代码有点错误，我现在是想用最小二乘法拟合自由表面的曲线，首先轮循所有自由表面网格，计算自由表面网格数量N，开辟数值存储自由表面上网格的重心（优良拟合三次曲线方程），gt,GS,max三个函数是用来求三次曲线方程y=a+bx+cx2+dx*x*x四个系数。然后把求出的四个系数保存到C_UDMI中，用于与离散的激光光线求交。代码编译之后有错误，可能思路上都已经与udf编程思想向偏离了，求各路大神指点，重谢！！！！！！！！！！！至少400金币！！！！呵呵
#include "udf.h"

#include "sg.h"

#include "sg_mphase.h"

#include "flow.h"

#include <math.h>

double gt(double x[],double y[]);
double GS(double a[4][4],double b[4],double x[4]);
double max(double array[4]);

DEFINE_ADJUST(store_gradient,domain)
{
  int phase_domain_index = 0;
  real xc[ND_ND],xd[ND_ND],x[4],n1,n2,thita_l,thita_g;
  Thread *t;
  Thread **pt;
  cell_t c;

  Domain *pDomain = DOMAIN_SUB_DOMAIN(domain,phase_domain_index);
  {
Alloc_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
Scalar_Reconstruction(pDomain, SV_VOF,-1,SV_VOF_RG,NULL);
Scalar_Derivatives(pDomain,SV_VOF,-1,SV_VOF_G,SV_VOF_RG,Vof_Deriv_Accumulate);
  }
  mp_thread_loop_c(t,domain,pt)
{
  if (FLUID_THREAD_P(t))
  {
int N=0;
Thread *ppt = pt[phase_domain_index];

begin_c_loop (c,t)
{
   C_CENTROID(xd,c,t);
   if(C_VOF(c, ppt)>0.1&&C_VOF(c, ppt)<0.9&&fabs(xd[0]-0)>0.0005)
      N++;
}
end_c_loop (c,t)

double a[N],b[N],c[N],d[N],e[N],f[N],g[N];

int i;
for(i=0; i<N; i++)
{
   a=1;
}

int j=0;
begin_c_loop (c,t)
{
   C_CENTROID(xc,c,t);
   if(C_VOF(c, ppt)>0.1&&C_VOF(c, ppt)<0.9&&fabs(xc[0]-0)>0.0005)
   {
      b[j]=xc[0];
      e[j]=xc[1];
      c[j]=b[j]*b[j];
      d[j]=b[j]*b[j]*b[j];
j++;
   }
}
end_c_loop (c,t)
f[0]=gt(a,a);f[1]=gt(a,b);f[2]=gt(a,c);f[3]=gt(a,d);f[4]=gt(b,b);f[5]=gt(b,c);f[6]=gt(b,d);
f[7]=gt(c,c);f[8]=gt(c,d);f[9]=gt(d,d);g[0]=gt(a,e);g[1]=gt(b,e);g[2]=gt(c,e);g[3]=gt(d,e);
double array[4][4]={{f[0],f[1],f[2],f[3]},{f[1],f[4],f[5],f[6]},{f[2],f[5],f[7],f[8]},{f[3],f[6],f[8],f[9]}};
GS(array,g,x);
C_UDMI(c,t,0) = x[0];
C_UDMI(c,t,1) = x[1];
C_UDMI(c,t,2) = x[2];
C_UDMI(c,t,3) = x[3];
  }
}
  Free_Storage_Vars(pDomain,SV_VOF_RG,SV_VOF_G,SV_NULL);
}

double gt(double x[],double y[])
{
double sum=0;
int i;
for(i=0;i<N;i++)
   sum=x*y+sum;
return sum;
}

double GS(double a[4][4],double b[4],double x[4])
{
double c[4]={0};
double x0[4]={0};
int i,j,k;
double sum=0;
for(k=1;;k++)
{
      for(i=0;i<4;i++)
      {
         for(j=0;j<4;j++)
         {
            sum=a[j]*x0[j]+sum;
         }
         x=x0+(b-sum)/a;
      c=fabs(x-x0);
      x0=x;
      sum=0;
}
r=max(c);
if(r<0.0001)
{
         break;
   }
}
return 0;
}

double max(double array[4])
{
double a=array[0];
int i;
for(i=0;i<4;i++)
{
if(a<array)
   a=array;
}
return a;
}