stream.c

来自「FreeFem++可以生成高质量的有限元网格。可以用于流体力学」· C语言 代码 · 共 1,962 行 · 第 1/4 页

  cb[2] = aire3 * dd;

  pt->cpt++;
  return(1);
}


/* return size of tetra */
double sizeTetra(pMesh mesh,int k) {
  pTetra   pt;
  pPoint   p[4];
  double   ax,ay,az,dd;
  double   hmin;
  int      i;
  static int  idire[6][2] = {{0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3}};

  pt = &mesh->tetra[k];
  for (i=0; i<4; i++)
    p[i] = &mesh->point[pt->v[i]];
  hmin = FLT_MAX;
  for (i=1; i<6; i++) {
    ax = p[idire[i][0]]->c[0] - p[idire[i][1]]->c[0];
    ay = p[idire[i][0]]->c[1] - p[idire[i][1]]->c[1];
    az = p[idire[i][0]]->c[2] - p[idire[i][1]]->c[2];
    dd = ax*ax + ay*ay + az*az;
    hmin = min(dd,hmin);
  }
  return(sqrt(hmin));
}

double sizeHexa(pMesh mesh,int k) {
  pHexa    ph;
  pPoint   p[8];
  double   ax,ay,az,dd;
  double   hmin;
  int      i;
  static int idire[12][2] = {{0,1}, {1,2}, {2,3}, {0,3}, {4,5}, {5,6},
                             {6,7}, {4,7}, {0,4}, {1,5}, {2,6}, {3,7}};

  ph = &mesh->hexa[k];
  for (i=0; i<8; i++)
    p[i] = &mesh->point[ph->v[i]];
  hmin = FLT_MAX;
  for (i=1; i<12; i++) {
    ax = p[idire[i][0]]->c[0] - p[idire[i][1]]->c[0];
    ay = p[idire[i][0]]->c[1] - p[idire[i][1]]->c[1];
    az = p[idire[i][0]]->c[2] - p[idire[i][1]]->c[2];
    dd = ax*ax + ay*ay + az*az;
    hmin = min(dd,hmin);
  }
  return(sqrt(hmin));
}

double sizeTria(pMesh mesh,int k) {
  pTriangle pt;
  pPoint    p0,p1;
  double    ax,ay,dd;
  double    hmin;
  int       i;
  static int idir[5] = {0,1,2,0,1};

  pt = &mesh->tria[k];
  hmin = FLT_MAX;
  for (i=0; i<3; i++) {
    p0 = &mesh->point[pt->v[i]];
    p1 = &mesh->point[pt->v[idir[i+1]]];
    ax = p0->c[0] - p1->c[0];
    ay = p0->c[1] - p1->c[1];
    dd = ax*ax + ay*ay;
    hmin = min(dd,hmin);
  }
  return(sqrt(hmin));
}

double sizeQuad(pMesh mesh,int k) {
  pQuad     pq;
  pPoint    p0,p1;
  double    ax,ay,dd;
  double    hmin;
  int       i;
  static int idir[7] = {0,1,2,3,0,1,2};

  pq = &mesh->quad[k];
  hmin = FLT_MAX;
  for (i=0; i<4; i++) {
    p0 = &mesh->point[pq->v[i]];
    p1 = &mesh->point[pq->v[idir[i+1]]];
    ax = p0->c[0] - p1->c[0];
    ay = p0->c[1] - p1->c[1];
    dd = ax*ax + ay*ay;
    hmin = min(dd,hmin);
  }
  return(sqrt(hmin));
}

/* vector interpolation */
double field3DInterp(pMesh mesh,int iel,double *cb,double *v) {
  pTetra     pt;
  pSolution  ps0,ps1,ps2,ps3;
  double     dd;

  pt  = &mesh->tetra[iel];
  ps0 = &mesh->sol[pt->v[0]];
  ps1 = &mesh->sol[pt->v[1]];
  ps2 = &mesh->sol[pt->v[2]];
  ps3 = &mesh->sol[pt->v[3]];

  v[0] = cb[0]*ps0->m[0] + cb[1]*ps1->m[0] + \
         cb[2]*ps2->m[0] + cb[3]*ps3->m[0];
  v[1] = cb[0]*ps0->m[1] + cb[1]*ps1->m[1] + \
         cb[2]*ps2->m[1] + cb[3]*ps3->m[1];
  v[2] = cb[0]*ps0->m[2] + cb[1]*ps1->m[2] + \
         cb[2]*ps2->m[2] + cb[3]*ps3->m[2];
  dd = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
  if ( dd > 0.0f ) {
    v[0] /= dd;
    v[1] /= dd;
    v[2] /= dd;
  }
  return(dd);
}

/* vector interpolation */
double vector3DInterp(pMesh mesh,pPoint pt[4],double *cb,double *v) {
  pSolution  ps0,ps1,ps2,ps3;
  double     dd;

  ps0 = &mesh->sol[pt[0]-&mesh->point[0]];
  ps1 = &mesh->sol[pt[1]-&mesh->point[0]];
  ps2 = &mesh->sol[pt[2]-&mesh->point[0]];
  ps3 = &mesh->sol[pt[3]-&mesh->point[0]];

  v[0] = cb[0]*ps0->m[0] + cb[1]*ps1->m[0] + \
         cb[2]*ps2->m[0] + cb[3]*ps3->m[0];
  v[1] = cb[0]*ps0->m[1] + cb[1]*ps1->m[1] + \
         cb[2]*ps2->m[1] + cb[3]*ps3->m[1];
  v[2] = cb[0]*ps0->m[2] + cb[1]*ps1->m[2] + \
         cb[2]*ps2->m[2] + cb[3]*ps3->m[2];
  dd = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
  if ( dd > 0.0f ) {
    v[0] /= dd;
    v[1] /= dd;
    v[2] /= dd;
  }
  return(dd);
}


double field2DInterp(pMesh mesh,int iel,double *cb,double *v) {
  pTriangle  pt;
  pSolution  ps0,ps1,ps2;
  double     dd;

  pt  = &mesh->tria[iel];
  ps0 = &mesh->sol[pt->v[0]];
  ps1 = &mesh->sol[pt->v[1]];
  ps2 = &mesh->sol[pt->v[2]];

  v[0] = cb[0]*ps0->m[0] + cb[1]*ps1->m[0] + cb[2]*ps2->m[0];
  v[1] = cb[0]*ps0->m[1] + cb[1]*ps1->m[1] + cb[2]*ps2->m[1];
  v[2] = 0.0f;
  dd = sqrt(v[0]*v[0] + v[1]*v[1]);
  if ( dd > 0.0f ) {
    v[0] /= dd;
    v[1] /= dd;
  }
  return(dd);
}


/* add point to display list, if needed */
int filterPoint(pScene sc,Stream *st,float *p,ubyte color) {
  double  norm,rgb[3],kc,ux,uy,uz,vx,vy,vz,dd;
  int     i;
  static double hsv[3] = { 0.0f, 1.0f, 0.80f };

  /* store point */
  memcpy(st->stpt[++st->stnp],p,3*sizeof(float));
  nbar++;

  /* point color */
  norm = st->norm;
  if ( !color ) {
    if ( norm < sc->iso.val[0] ) 
      norm = sc->iso.val[0];  
    else if ( norm > sc->iso.val[MAXISO-1] )
      norm = sc->iso.val[MAXISO-1];
    for (i=0; i<MAXISO-1; i++)
      if ( norm < sc->iso.val[i] )  break;
    kc = (norm-sc->iso.val[i-1]) / (sc->iso.val[i] - sc->iso.val[i-1]);
    st->stcol[st->stnp] = sc->iso.col[i-1]*(1.0-kc)+sc->iso.col[i]*kc;
    st->stiso[st->stnp] = i;
  }
  else {
    st->stcol[st->stnp] = sc->iso.col[MAXISO-1];
    st->stiso[st->stnp] = MAXISO-1;
  }

  if ( sc->mode & S_ALTITUDE )
    st->stpt[st->stnp][2] = altcoef*norm;

  if ( !color && (st->stiso[0] != st->stiso[st->stnp]) ) {
    hsv[0] = st->stcol[st->stnp];
    hsvrgb(hsv,rgb);
    glColor3dv(rgb);
    glVertex3fv(st->stpt[st->stnp]);

    memcpy(st->stpt[0],st->stpt[st->stnp],3*sizeof(float));
    st->stcol[0] = st->stcol[st->stnp];
    st->stiso[0] = st->stiso[st->stnp];
    st->stnp = 0;
    return(1);
  }

  if ( st->stnp < 2 )  return(0);

  /* filtering point */
  ux = st->stpt[0][0] - st->stpt[1][0];
  uy = st->stpt[0][1] - st->stpt[1][1];
  uz = st->stpt[0][2] - st->stpt[1][2];
  dd = ux*ux + uy*uy + uz*uz;
  
  if ( dd > 0.0 ) {  
    dd = 1.0f / sqrt(dd);
    ux *= dd;
    uy *= dd;
    uz *= dd;
  }
  else {
    memcpy(st->stpt[0],st->stpt[1],2*3*sizeof(float));
    st->stnp--;
    return(0);
  }

  vx = st->stpt[2][0] - st->stpt[1][0];
  vy = st->stpt[2][1] - st->stpt[1][1];
  vz = st->stpt[2][2] - st->stpt[1][2];
  dd = vx*vx + vy*vy + vz*vz;
  if ( dd > 0.0 ) {
    dd = 1.0f / sqrt(dd);
    vx *= dd;
    vy *= dd;
    vz *= dd;
  }
  else {
    memcpy(st->stpt[1],st->stpt[2],3*sizeof(float));
    st->stnp--;
    return(0);
  }

  dd = ux*vx + uy*vy + uz*vz;
  if ( dd > COS178 ) {
    hsv[0] = st->stcol[st->stnp];
    hsvrgb(hsv,rgb);
    glColor3dv(rgb);
    glVertex3fv(st->stpt[st->stnp]);

    memcpy(st->stpt[0],st->stpt[st->stnp],3*sizeof(float));
    st->stcol[0] = st->stcol[st->stnp];
    st->stiso[0] = st->stiso[st->stnp];
    st->stnp = 0;
    return(1);
  }
  else {
    memcpy(st->stpt[st->stnp-1],st->stpt[st->stnp],3*sizeof(float));
    st->stcol[st->stnp-1] = st->stcol[st->stnp];
    st->stiso[st->stnp-1] = st->stiso[st->stnp];
    st->stnp--;
    return(0);
  }
}


/* add vertex to display list */
void addPoint(pScene sc,Stream *st,float *p,ubyte color) {
  double  norm,rgb[3],kc;
  int     i;
  static double hsv[3] = { 0.0f, 1.0f, 0.80f };

  /* point color */
  norm = st->norm;
  i = MAXISO-1;
  if ( !color ) {
    norm = st->norm;
    if ( norm < sc->iso.val[0] )
      norm = sc->iso.val[0];  
    else if ( norm > sc->iso.val[MAXISO-1] )
      norm = sc->iso.val[MAXISO-1];
    for (i=0; i<MAXISO-1; i++)
      if ( norm < sc->iso.val[i] )  break;
    kc = (norm-sc->iso.val[i-1]) / (sc->iso.val[i] - sc->iso.val[i-1]);
    hsv[0] = sc->iso.col[i-1]*(1.0-kc)+sc->iso.col[i]*kc;
  }

  if ( sc->mode & S_ALTITUDE )
    p[2] = altcoef*norm;

  hsvrgb(hsv,rgb);
  glColor3dv(rgb);
  glVertex3fv(p);
  st->stnp = 0;
  memcpy(st->stpt[st->stnp],p,3*sizeof(float));
  st->stcol[st->stnp] = hsv[0];
  st->stiso[st->stnp] = i;
  st->stnp++;
}


int nxtPoint3D(pMesh mesh,int nsdep,float *p,float step,double *v) {
  pTetra     pt;
  double     norm,h6,cb[4],v1[3],v2[3],v3[3];
  float      xp1[3],xp2[3],xp3[3];
  int        k;

  /* 4th order Runge-Kutta */
  xp1[0] = p[0] + 0.5*step*v[0];
  xp1[1] = p[1] + 0.5*step*v[1];
  xp1[2] = p[2] + 0.5*step*v[2];

  k = locateTetra(mesh,nsdep,++mesh->mark,xp1,cb);
  if ( !k )  return(0);
  norm = field3DInterp(mesh,k,cb,v1);
  pt = &mesh->tetra[k];
  pt->cpt--;

  xp2[0] = p[0] + 0.5*step*v1[0];
  xp2[1] = p[1] + 0.5*step*v1[1];
  xp2[2] = p[2] + 0.5*step*v1[2];

  k = locateTetra(mesh,k,++mesh->mark,xp2,cb);
  if ( !k )  return(0);
  norm = field3DInterp(mesh,k,cb,v2);
  pt = &mesh->tetra[k];
  pt->cpt--;

  xp3[0] = p[0] + step*v2[0];
  xp3[1] = p[1] + step*v2[1];
  xp3[2] = p[2] + step*v2[2];

  k = locateTetra(mesh,k,++mesh->mark,xp3,cb);
  if ( !k )  return(0);
  norm = field3DInterp(mesh,k,cb,v3);
  pt = &mesh->tetra[k];
  pt->cpt--;

  h6    = step / 6.0;
  p[0] += h6 * (v[0] + 2*(v1[0] + v2[0]) + v3[0]);
  p[1] += h6 * (v[1] + 2*(v1[1] + v2[1]) + v3[1]);
  p[2] += h6 * (v[2] + 2*(v1[2] + v2[2]) + v3[2]);

  return(1);
}

int nxtPoint2D(pMesh mesh,int nsdep,float *p,float step,double *v) {
  pTriangle  pt;
  double     norm,h6,cb[3],v1[3],v2[3],v3[3];
  float      xp1[3],xp2[3],xp3[3];
  int        k;

  /* 4th order Runge-Kutta */
  xp1[0] = p[0] + 0.5*step*v[0];
  xp1[1] = p[1] + 0.5*step*v[1];

  k = locateTria(mesh,nsdep,++mesh->mark,xp1,cb);
  if ( !k )  return(0);
  norm = field2DInterp(mesh,k,cb,v1);
  pt = &mesh->tria[k];
  pt->cpt--;
  
  xp2[0] = p[0] + 0.5*step*v1[0];
  xp2[1] = p[1] + 0.5*step*v1[1];

  k = locateTria(mesh,k,++mesh->mark,xp2,cb);
  if ( !k )  return(0);
  norm = field2DInterp(mesh,k,cb,v2);
  pt = &mesh->tria[k];
  pt->cpt--;

  xp3[0] = p[0] + step*v2[0];
  xp3[1] = p[1] + step*v2[1];

  k = locateTria(mesh,k,++mesh->mark,xp3,cb);
  if ( !k )  return(0);
  norm = field2DInterp(mesh,k,cb,v3);
  pt = &mesh->tria[k];
  pt->cpt--;

  h6 = step / 6.0;
  p[0] += h6 * (v[0] + 2*(v1[0] + v2[0]) + v3[0]);
  p[1] += h6 * (v[1] + 2*(v1[1] + v2[1]) + v3[1]);

  return(1);
}


/* read streamlines origins */
int parseStream(pScene sc,pMesh mesh) {
  FILE    *in;
  pStream  st;
  float    x,y,z;
  int      i,k,nbp,ret;
  char    *ptr,data[128],key[256],tmp[128];

  /* input file */
  strcpy(tmp,mesh->name);
  ptr = (char*)strstr(tmp,".mesh");
  if ( ptr )  *ptr = '\0';

  sprintf(data,"%s.iso",tmp);
  in = fopen(data,"r");
  if ( !in ) {
    sscanf(data,"DEFAULT.iso");
    in = fopen(data,"r");
    if ( !in )  return(0);
  }
  if ( !quiet )  fprintf(stdout,"  Reading %s\n",data);
  sc->stream = createStream(sc,mesh); 
  st = sc->stream;

  while ( !feof(in) ) {
    fscanf(in,"%s",key);
    for (i=0; i<strlen(key); i++) key[i] = tolower(key[i]);

    if ( !strcmp(key,"nblines") ) {
      fscanf(in,"%d",&nbp);
      st->nbstl = nbp;
      if ( mesh->dim == 3 )
        for (k=1; k<=3*st->nbstl; k+=3) {
          ret = fscanf(in,"%f %f %f\n",&x,&y,&z);
printf("x %f %f %f\n",x,y,z);
          st->listp[k]   = x - mesh->xtra;
          st->listp[k+1] = y - mesh->ytra;
          st->listp[k+2] = z - mesh->ztra;
printf("x %f %f %f\n",st->listp[k],st->listp[k+1],st->listp[k+2]);
        }
      else
        for (k=1; k<=2*st->nbstl; k+=2) {
          ret = fscanf(in,"%f %f\n",&x,&y);
          st->listp[k]   = x - mesh->xtra;
          st->listp[k+1] = y - mesh->ytra;
        }
    }
    else if ( !strcmp(key,"euler") ) {
      st->typtrack = Euler;
    }
    else if ( !strcmp(key,"box") ) {
      ret = fscanf(in,"%f %f",&x,&y);
      if ( ret != 2 )  break; 
      st->xmin = x - mesh->xtra;
      st->xmax = y - mesh->xtra;
      ret = fscanf(in,"%f %f",&x,&y);
      if ( ret != 2 )  break;
      st->ymin = x - mesh->ytra;
      st->ymax = y - mesh->ytra; 
      if ( mesh->dim == 3 ) {
        ret = fscanf(in,"%f %f",&x,&y);
        if ( ret != 2 )  break;
        st->zmin = x - mesh->ztra;
        st->zmax = y - mesh->ztra;
      }
    }
    else if ( key[0] == '#' ) {
      fgets(key,255,in);
    }
  }
  fclose(in);

  if ( !st->nbstl ) {
    fprintf(stderr,"   ## No data found.\n");
    return(0);
  }
k = 1;
printf("fin proc %f %f %f\n",sc->stream->listp[k],sc->stream->listp[k+1],sc->stream->listp[k+2]);

  return(1);
}


/* build lists for streamlines */
int listTetraStream(pScene sc,pMesh mesh,float *pp,int squiet) {
  pTetra     pt;
  pStream    st;
  double     dd,cb[4],cbdep[4],v[4],vdep[4],sizedep,normdep;
  float      step,p[3],ox,oy,oz,ldt;
  int        i,k,exh,depart,nsdep,nsfin,nsold,nbp,maxpts;
  clock_t    ct;
  FILE      *out;

  /* default */
  if ( !mesh->ntet )  return(0);
  else if ( egal(sc->iso.val[0],sc->iso.val[MAXISO-1]) )  return(0);