stream.c

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

  if ( ddebug ) printf("\n create streamlines list / TETRA\n");

  if ( !squiet && !ddebug ) {
    fprintf(stdout," Building streamline(s)");
    fflush(stdout);
  }
  ct = clock();

  /* build display list */
  st = sc->stream;
  if ( st->nbstl > MAX_LST-1 )  return(0);
  sc->slist[st->nbstl] = glGenLists(1);
  glNewList(sc->slist[st->nbstl],GL_COMPILE);
  if ( glGetError() )  return(0);

  k = st->nbstl*3 + 1;
  st->listp[k+0] = pp[0];
  st->listp[k+1] = pp[1];
  st->listp[k+2] = pp[2];
  st->nbstl++;
printf("\n%d: pp = %f %f %f\n",st->nbstl,st->listp[k+0],st->listp[k+1],st->listp[k+2]);

  maxpts = max(MAX_PTS,5*mesh->ntet);
  glLineWidth(2.0);

  /* compute streamline */
  nbp   = 0;
  exh   = 0;
  nsdep = mesh->ntet / 2;
  step  = 0.0;
  nbar  = 0;
  if ( ddebug ) 
    printf("   start point %d: %f %f %f\n",
           3*k/3,st->listp[k],st->listp[k+1],st->listp[k+2]);

  for (i=1; i<mesh->ntet; i++) {
    pt = &mesh->tetra[i];
    pt->cpt = 0;
  }

  /* find enclosing tet */
  memcpy(p,pp,3*sizeof(float));
  depart = locateTetra(mesh,nsdep,++mesh->mark,p,cb);
printf("depart = %d\n",depart);
  if ( !depart ) {
    for (depart=1; depart<=mesh->ntet; depart++) {
      pt = &mesh->tetra[depart];
      if ( pt->mark != mesh->mark && inTetra(mesh,depart,p,cb) )
        break;
    }
    if ( depart > mesh->ntet ) {
      glEndList();
      return(0);
    }
  }

  st->norm = field3DInterp(mesh,depart,cb,v);
  memcpy(cbdep,cb,4*sizeof(double));
  memcpy(vdep,v,4*sizeof(double));
  st->size = sizeTetra(mesh,depart);
  sizedep  = st->size;
  normdep  = st->norm;
  ldt      = 0.0;

  if ( st->size == 0.0 )
    step = EPS*sc->dmax;
  else
    step = HSIZ * min(st->size,st->norm);

  /* build display list incrementally */
  nsdep = nsold = depart;
  glBegin(GL_LINE_STRIP);
  addPoint(sc,st,p,0);
  nbp++;

  if ( sc->par.maxtime < FLT_MAX ) {
    sc->par.cumtim = 0.0;
    step = min(0.05*sc->par.dt,step);
    out  = fopen("particules.dat","a+");
    assert(out);
    fprintf(out,"\n%8.2f  %f %f %f\n",
            sc->par.cumtim,p[0]+mesh->xtra,p[1]+mesh->ytra,p[2]+mesh->ztra);
  }

  do {
    ox = p[0];
    oy = p[1];
    oz = p[2];

    /* move to next point */
    if ( st->typtrack == Euler || !nxtPoint3D(mesh,nsdep,p,step,v) ) {
      p[0] += step*v[0];
      p[1] += step*v[1];
      p[2] += step*v[2];
    }
    if ( p[0]<st->xmin || p[1]<st->ymin || p[2]<st->zmin ||
         p[0]>st->xmax || p[1]>st->ymax || p[2]>st->zmax )
      break;
    else if ( sc->par.maxtime < FLT_MAX ) {
      ox -= p[0];
      oy -= p[1];
      oz -= p[2];
      dd  = sqrt(ox*ox + oy*oy + oz*oz) / st->norm;
      ldt += dd;
      sc->par.cumtim     += dd;
      if ( sc->par.cumtim  > sc->par.maxtime )  break;
      
      if ( ldt > sc->par.dt  ) {
        fprintf(out,"%8.2f  %f %f %f\n",
	        sc->par.cumtim,p[0]+mesh->xtra,p[1]+mesh->ytra,p[2]+mesh->ztra);
    	ldt = fabs(sc->par.dt - ldt);
      }
    }

    /* find tet containing p */
    nsfin = locateTetra(mesh,nsdep,++mesh->mark,p,cb);
    if ( !nsfin )  break;
    nsdep = nsfin;
    pt = &mesh->tetra[nsdep];
    if ( pt->cpt > MAX_CPT )  break;

    /* adjust local stepsize */
    if ( nsdep != nsold ) {
      st->size = sizeTetra(mesh,nsdep);
      nsold = nsdep;
    }

    /* vector field interpolation */
    st->norm = field3DInterp(mesh,nsdep,cb,v);
    step     = HSIZ*min(st->size,st->norm);
    if ( sc->par.maxtime < FLT_MAX )
      step = min(0.05*sc->par.dt,step);
    if ( step == 0.0 )  break;

    nbp += filterPoint(sc,st,p,0);
  }
  while ( nbp < maxpts );
  addPoint(sc,st,p,0);
  glEnd();

  if ( nbp >= maxpts || sc->par.maxtime < FLT_MAX ) {
    glLineWidth(1.0);
    glEndList();
    if ( !squiet && !ddebug ) {
      fprintf(stdout,": %d (%d, %.2f) / %d lines",
              nbar,nbp,(float)nbp/nbar,k/3);
      ct = difftime(clock(),ct);
      fprintf(stdout," %6.2f sec.\n",ct);
    }
    if ( sc->par.maxtime < FLT_MAX ) {
      fprintf(out,"%8.2f  %f %f %f\n",
              sc->par.cumtim,p[0]+mesh->xtra,p[1]+mesh->ytra,p[2]+mesh->ztra);
      fclose(out);
    }
    return(1);
  }

  /* reverse orientation */
  memcpy(p,pp,3*sizeof(float));
  memcpy(cb,cbdep,4*sizeof(double));
  memcpy(v,vdep,4*sizeof(double));
  st->norm = normdep;
  st->size = sizedep;
  if ( st->size == 0.0 )
    step = EPS * sc->dmax;
  else
    step = HSIZ * min(st->size,st->norm);

  /* build display list incrementally */
  nsdep = nsold = depart;
  glBegin(GL_LINE_STRIP);
  addPoint(sc,st,p,0);
  nbp++;

  do {
    /* move to next point */
    if ( st->typtrack == Euler || !nxtPoint3D(mesh,nsdep,p,-step,v) ) {
      p[0] -= step*v[0];
      p[1] -= step*v[1];
      p[2] -= step*v[2];
    }
    if ( p[0]<st->xmin || p[1]<st->ymin || p[2]<st->zmin ||
         p[0]>st->xmax || p[1]>st->ymax || p[2]>st->zmax )
      break;

    /* find tet containing p */
    nsfin = locateTetra(mesh,nsdep,++mesh->mark,p,cb);
    if ( !nsfin )  break;
    nsdep = nsfin;
    pt = &mesh->tetra[nsdep];
    if ( pt->cpt > MAX_CPT )  break;

    /* adjust local stepsize */
    if ( nsdep != nsold ) {
      st->size = sizeTetra(mesh,nsdep);
      nsold = nsdep;
    }

    /* vector field interpolation */
    st->norm = field3DInterp(mesh,nsdep,cb,v);
    step = HSIZ * min(st->size,st->norm);
    if ( step == 0.0 )  break;

    nbp += filterPoint(sc,st,p,0);
  }
  while ( nbp < maxpts );
  addPoint(sc,st,p,0);
  glEnd();
  glLineWidth(1.0);
  glEndList();

  if ( !nbp ) {
    st->nbstl--;
    if ( !squiet && !ddebug )  fprintf(stdout,"..empty\n");
    return(0);
  }

  if ( !squiet && !ddebug ) {  
    if ( nbar )
      fprintf(stdout,": %d (%d, %.2f) / %d lines",nbar,nbp,(float)nbp/nbar,k/3);
    ct = difftime(clock(),ct);
    fprintf(stdout," %6.2f sec.\n",ct);
  }

  return(1);
}


int listHexaStream(pScene sc,pMesh mesh,float *pp,int squiet) {
  pHexa      ph;
  pPoint     pt[4];
  pStream    st;
  double     cbdep[4],cb[4],v[6],vdep[6],sizedep,normdep;
  float      step,p[3];
  int        i,k,exh,depart,nsdep,nsfin,nsold,nbp,maxpts;
  mytime     tt;

  /* default */
  if ( !mesh->nhex )  return(0);
  else if ( egal(sc->iso.val[0],sc->iso.val[MAXISO-1]) )  return(0);
  if ( ddebug ) printf("create streamlines list / HEXA\n");

  if ( !squiet && !ddebug ) {
    fprintf(stdout," Building streamline(s)");
    fflush(stdout);
    chrono(ON,&tt);
  }

  /* build display list */
  st = sc->stream;
  if ( st->nbstl > MAX_LST-1 )  return(0);
  sc->slist[st->nbstl] = glGenLists(1);
  glNewList(sc->slist[st->nbstl],GL_COMPILE);
  if ( glGetError() )  return(0);

  st->nbstl++;
  k = st->nbstl*3;
  st->listp[k]   = pp[0];
  st->listp[k+1] = pp[1];
  st->listp[k+2] = pp[2];

  maxpts = max(MAX_PTS,5*mesh->nhex);
  glLineWidth(2.0);

  /* compute streamline */
  nbp   = 0;
  exh   = 0;
  nsdep = mesh->nhex / 2;
  step  = 0.0f;
  nbar  = 0;
  if ( ddebug ) 
    printf("   start point %d: %f %f %f\n",
           3*k/3,st->listp[k],st->listp[k+1],st->listp[k+2]);

  for (i=1; i<mesh->nhex; i++) {
    ph = &mesh->hexa[i];
    ph->cpt  = 0;
  }

  /* find enclosing tet */
  memcpy(p,pp,3*sizeof(float));
  depart = locateHexa(mesh,nsdep,++mesh->mark,p,cb,pt);
printf("DEPART %d\n",depart);
ph = &mesh->hexa[depart];
printf("sommets %d %d %d %d %d %d %d %d\n",
ph->v[0],ph->v[1],ph->v[2],ph->v[3],ph->v[4],ph->v[5],ph->v[6],ph->v[7]);
  
  if ( !depart ) {
    for (depart=1; depart<=mesh->nhex; depart++) {
      ph = &mesh->hexa[depart];
      if ( ph->mark != mesh->mark && inHexa(mesh,depart,p,cb,pt) )
        break;
    }
    if ( depart > mesh->nhex )  return(0);
  }

  st->norm = vector3DInterp(mesh,pt,cb,v);
  memcpy(cbdep,cb,4*sizeof(double));
  memcpy(vdep,v,4*sizeof(double));
  st->size = sizeHexa(mesh,depart);
  sizedep  = st->size;
  normdep  = st->norm;
  if ( st->size == 0.0f )
    step = EPS*sc->dmax;
  else
    step = HSIZ * st->size;

  /* build display list incrementally */
  nsdep = nsold = depart;
  glBegin(GL_LINE_STRIP);
  addPoint(sc,st,p,0);
  nbp++;

  do {
    /* move to next point */
    if ( st->typtrack == Euler || !nxtPoint3D(mesh,nsdep,p,step,v)) {
      p[0] += step*v[0];
      p[1] += step*v[1];
      p[2] += step*v[2];
    }
    if ( p[0]<st->xmin || p[1]<st->ymin || p[2]<st->zmin ||
         p[0]>st->xmax || p[1]>st->ymax || p[2]>st->zmax )
      break;

    /* find tet containing p */
    nsfin = locateHexa(mesh,nsdep,++mesh->mark,p,cb,pt);
    if ( !nsfin )  break;
    nsdep = nsfin;
    ph = &mesh->hexa[nsdep];
    if ( ph->cpt > MAX_CPT )  break;

    /* adjust local stepsize */
    if ( nsdep != nsold ) {
      st->size = sizeHexa(mesh,nsdep);
      step  = HSIZ*st->size;
      nsold = nsdep;
    }

    /* vector field interpolation */
    st->norm = vector3DInterp(mesh,pt,cb,v);
    if ( st->norm < EPS*step )  break;
    step = min(step,st->norm);
    if ( step == 0.0f )  break; /*step = 1.0e-06*sc->dmax;*/

    nbp += filterPoint(sc,st,p,0);
  }
  while ( nbp < maxpts );
  glEnd();

  if ( nbp >= maxpts ) {
    glLineWidth(1.0);
    glEndList();
    if ( !squiet && !ddebug ) {
      fprintf(stdout,": %d (%d, %.2f) / %d lines",
              nbar,nbp,(float)nbp/nbar,k/3);
      chrono(OFF,&tt);
      fprintf(stdout," %6.2f sec.\n",gttime(tt));
    }
    return(1);
  }

  /* reverse orientation */
  memcpy(p,pp,3*sizeof(float));
  memcpy(cb,cbdep,4*sizeof(double));
  memcpy(v,vdep,4*sizeof(double));
  st->norm = normdep;
  st->size = sizedep;
  if ( st->size == 0.0f )
    step = EPS * sc->dmax;
  else
    step = HSIZ*st->size;

  /* build display list incrementally */
  nsdep = nsold = depart;
  glBegin(GL_LINE_STRIP);
  addPoint(sc,st,p,0);
  nbp++;

  do {
    /* move to next point */
    if ( st->typtrack == Euler || !nxtPoint3D(mesh,nsdep,p,-step,v) ) {
      p[0] -= step*v[0];
      p[1] -= step*v[1];
      p[2] -= step*v[2];
    }
    if ( p[0]<st->xmin || p[1]<st->ymin || p[2]<st->zmin ||
         p[0]>st->xmax || p[1]>st->ymax || p[2]>st->zmax )
      break;

    /* find tet containing p */
    nsfin = locateHexa(mesh,nsdep,++mesh->mark,p,cb,pt);
    if ( !nsfin )  break;
    nsdep = nsfin;
    ph = &mesh->hexa[nsdep];
    if ( ph->cpt > MAX_CPT )  break;

    /* adjust local stepsize */
    if ( nsdep != nsold ) {
      st->size = sizeHexa(mesh,nsdep);
      step  = HSIZ * st->size;
      nsold = nsdep;
    }

    /* vector field interpolation */
    st->norm = vector3DInterp(mesh,pt,cb,v);
    if ( st->norm < EPS*step )   break;
    step = min(step,st->norm);
    if ( step == 0.0f )  break; /*step = 1.e-06 * sc->dmax;*/

    nbp += filterPoint(sc,st,p,0);
  }
  while ( nbp < maxpts );
  glEnd();
  glLineWidth(1.0);
  glEndList();

  if ( nbp >= maxpts ) {
    glLineWidth(1.0);
    glEndList();
    if ( !squiet && !ddebug ) {
      fprintf(stdout,": %d (%d, %.2f) / %d lines",
              nbar,nbp,(float)nbp/nbar,k/3);
      chrono(OFF,&tt);
      fprintf(stdout," %6.2f sec.\n",gttime(tt));
    }
    return(1);
  }

  if ( !nbp ) {
    st->nbstl--;
    if ( !squiet && !ddebug )  fprintf(stdout,"..empty\n");
    return(0);
  }

  if ( !squiet && !ddebug ) {  
    if ( nbar )
      fprintf(stdout,": %d (%d, %.2f) / %d lines",nbar,nbp,(float)nbp/nbar,k/3);
    chrono(OFF,&tt);
    fprintf(stdout," %6.2f sec.\n",gttime(tt));
  }

  return(1);
}


int listTriaStream(pScene sc,pMesh mesh,float *pp) {
  pTriangle  pt;
  pStream    st;
  double     dd,cb[3],cbdep[3],v[3],vdep[3],sizedep,normdep;
  float      step,p[3],ox,oy,ldt;
  int        i,k,exh,depart,nsdep,nsfin,nsold,nbp,maxpts;
  clock_t    ct;
  FILE      *out;

  /* default */
  if ( !mesh->nt )  return(0);
  if ( ddebug ) printf("create streamlines / TRIA\n");
  if ( egal(sc->iso.val[0],sc->iso.val[MAXISO-1]) )  return(0);

  fprintf(stdout," Building streamline(s)");
  fflush(stdout);
  ct = clock();

  /* build display list */
  st = sc->stream;
  if ( st->nbstl > MAX_LST-1 )  return(0);
  sc->slist[st->nbstl] = glGenLists(1);
  glNewList(sc->slist[st->nbstl],GL_COMPILE);
  if ( glGetError() )  return(0);

  st->nbstl++;
  k = st->nbstl*3;
  st->listp[k]   = pp[0];
  st->listp[k+1] = pp[1];
  st->listp[k+2] = pp[2];

  maxpts = max(MAX_PTS,5*mesh->nt);
  glLineWidth(2.0);

  /* compute streamlines */
  nbp   = 0;
  exh   = 0;
  nsdep = mesh->nt / 2;
  step  = 0.0;
  nbar  = 0;
  if ( ddebug )
    printf("   start point %d: %f %f\n",3*k/3,st->listp[k],st->listp[k+1]);

  for (i=1; i<=mesh->nt; i++) {
    pt = &mesh->tria[i];
    pt->cpt = 0;