stream.c

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

  }

  /* find enclosing triangle */
  memcpy(p,pp,3*sizeof(float));
  depart = locateTria(mesh,nsdep,++mesh->mark,p,cb);
  if ( !depart ) {
    if ( ddebug )  printf("exhaustif search\n");
    for (depart=1; depart<=mesh->nt; depart++) {
      pt = &mesh->tria[depart];
      if ( pt->mark != mesh->mark && inTria(mesh,depart,p,cb) )
        break;
    }
    if ( depart > mesh->nt ) {
      st->nbstl--;
      glEndList();
      glLineWidth(1.0);
      fflush(stdout);
      return(0);
    }
  }

  st->norm = field2DInterp(mesh,depart,cb,v);
  memcpy(cbdep,cb,3*sizeof(double));
  memcpy(vdep,v,3*sizeof(double));
  st->size   = sizeTria(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\n",
            sc->par.cumtim,p[0]+mesh->xtra,p[1]+mesh->ytra);
  }

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

    /* move to next point */
    if ( st->typtrack == Euler || !nxtPoint2D(mesh,nsdep,p,step,v) ) {
      p[0] += step*v[0];
      p[1] += step*v[1];
    }
    if ( p[0]<st->xmin || p[1]<st->ymin ||
         p[0]>st->xmax || p[1]>st->ymax )
      break;
    else if ( sc->par.maxtime < FLT_MAX ) {
      ox -= p[0];
      oy -= p[1];
      dd  = sqrt(ox*ox + oy*oy) / 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\n",
                sc->par.cumtim,p[0]+mesh->xtra,p[1]+mesh->ytra);
        ldt = fabs(sc->par.dt - ldt);
      }
    }

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

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

    /* vector field interpolation */
    st->norm = field2DInterp(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();
    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\n",
              sc->par.cumtim,p[0]+mesh->xtra,p[1]+mesh->ytra);
      fclose(out);
    }
    return(1);
  }

  /* reverse orientation */
  memcpy(p,pp,3*sizeof(float));
  memcpy(cb,cbdep,3*sizeof(double));
  memcpy(v,vdep,3*sizeof(double));
  st->norm = normdep;
  st->size = sizedep;
  if ( st->size == 0.0 )
    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 || !nxtPoint2D(mesh,nsdep,p,-step,v) ) {
      p[0] -= step*v[0];
      p[1] -= step*v[1];
    }
    if ( p[0]<st->xmin || p[1]<st->ymin || p[0]>st->xmax || p[1]>st->ymax )
      break;

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

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

    /* vector field interpolation */
    st->norm = field2DInterp(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--;
    fprintf(stdout,".. empty.\n");
    return(0);
  }

  /*if ( 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 listSaddleStream(pScene sc,pMesh mesh,int depart,
                     float *pp,float *vv,double lambda) {
  pTriangle  pt;
  pStream    st;
  double     cb[3],v[3];
  float      sens,step,p[3];
  int        i,k,nsdep,nsfin,nsold,nbp,maxpts;

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

  /* 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);
  maxpts = max(MAX_PTS,5*mesh->nt);
  glLineWidth(2.0);

  st->nbstl++;
  k = st->nbstl*3;
  st->listp[k]   = p[0] = pp[0];
  st->listp[k+1] = p[1] = pp[1];
  st->listp[k+2] = p[2] = 0.0f;
  
  for (i=1; i<=mesh->nt; i++) {
    pt = &mesh->tria[i];
    pt->cpt = 0;
  }
 
  /* compute streamlines */
  nsold = nsdep = depart;
  nbp   = nbar   = 0;
  glBegin(GL_LINE_STRIP);
  addPoint(sc,st,pp,1);

  st->size = sizeTria(mesh,depart);
  st->norm = sqrt(vv[0]*vv[0] + vv[1]*vv[1]);
  if ( st->size == 0.0f )
    step = EPS * sc->dmax;
  else
    step = HSIZ* st->size;
  if ( st->norm > 0.0f ) {
    v[0] = vv[0] / st->norm;
    v[1] = vv[1] / st->norm;
  }
  sens = lambda < 0.0f ? -1. : 1.;

  /* build display list incrementally */
  do {
    /* move to next point */
    if ( st->typtrack == Euler || !nxtPoint2D(mesh,nsdep,p,step,v) ) {
      p[0] += sens*step*v[0];
      p[1] += sens*step*v[1];
    }

    if ( p[0]<st->xmin || p[1]<st->ymin ||
         p[0]>st->xmax || p[1]>st->ymax )  break;

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

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

    /* vector field interpolation */
    st->norm = field2DInterp(mesh,nsdep,cb,v);
    if ( st->norm < EPS*step )  break;
    step = min(step,st->norm);
    if ( step == 0.0f )  break;

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

  if ( !nbp ) {
    glDeleteLists(sc->slist[st->nbstl--],1);
    return(0);
  }

  return(1);
}


pStream createStream(pScene sc,pMesh mesh) {
  pStream    st;

  /* hash simplices */
  if ( mesh->dim == 2 ) {
    if ( mesh->nt && !hashTria(mesh) )    return(0);
  }
  else {
    if ( mesh->ntet && !hashTetra(mesh) ) return(0);
    if ( mesh->nhex && !hashHexa(mesh) )  return(0);
  }

  st = (pStream)calloc(1,sizeof(struct sstream));
  if ( !st )  return(0);

  /*st->typtrack = Euler;*/
  st->typtrack = RK4;
  st->stnp     = 0;
  st->nbstl    = 0;

  /* bounding box */
  st->xmin = mesh->xmin - mesh->xtra;
  st->ymin = mesh->ymin - mesh->ytra;
  st->zmin = mesh->zmin - mesh->ztra;
  st->xmax = mesh->xmax - mesh->xtra;
  st->ymax = mesh->ymax - mesh->ytra;
  st->zmax = mesh->zmax - mesh->ztra;

  /* init list */
  st->listp = (float*)malloc((MAX_LST*3+1)*sizeof(float));
  assert(st->listp);

  sc->slist = (GLuint*)calloc(MAX_LST,sizeof(GLuint));
  if ( !sc->slist )  return(0);

  return(st);
}

/* create from point picking */
int streamRefPoint(pScene sc,pMesh mesh) {
  pPoint    ppt;
  float     s[3];

  ppt = &mesh->point[refitem];
  if ( ppt->flag )  return(0);
  s[0] = ppt->c[0];
  s[1] = ppt->c[1];
  s[2] = ppt->c[2];
  if ( mesh->dim == 2 )
    listTriaStream(sc,mesh,s);
  else {
    if ( mesh->ntet )      listTetraStream(sc,mesh,s,0);
    else if ( mesh->nhex ) listHexaStream(sc,mesh,s,0);
  }
  ppt->flag = 1;

  return(1);
}


/* read starting point in file.iso */
int streamIsoPoint(pScene sc,pMesh mesh) {
  pStream   st;
  int       k,nbp,nbstl;
  time_t    t;

  if ( !parseStream(sc,mesh) )  return(0);
  t = clock();
  fprintf(stdout," Building streamline(s)");
  fflush(stdout);

  st    = sc->stream;
  nbstl = st->nbstl;
  nbp   = 0;
  st->nbstl = 0;
  if ( mesh->dim == 3 ) {
    if ( !mesh->ntet )  return(0);
    nbp = 0;
    for (k=1; k<=3*nbstl; k+=3) {
printf("\n ici: %f %f %f\n",st->listp[k],st->listp[k+1],st->listp[k+2]);
      nbp += listTetraStream(sc,mesh,&st->listp[k],1);
    }
  }

  if ( !nbp )    return(0);
  if ( ddebug )  printf("stream start: %d points  ",nbp);
  fprintf(stdout,": %d lines",nbp);
  t = clock() - t;
  fprintf(stdout," %6.2f sec.\n",t/(float)CLOCKS_PER_SEC);

  return(1);
}


int streamRefTria(pScene sc,pMesh mesh) {
  pMaterial   pm;
  pTriangle   pt;
  pPoint      ppt;
  float       s[3];
  int         i,k,nmat,nbp,base;
  time_t      t;

  /* build list */
  base = ++mesh->mark;
  nbp  = 0;
  pt   = &mesh->tria[refitem];

  nmat = matRef(sc,pt->ref);
  /*nmat = !pt->ref ? DEFAULT_MAT : 1+(pt->ref-1)%(sc->par.nbmat-1);*/
  pm = &sc->material[nmat];
  k  = pm->depmat[LTria];
  if ( !k || pm->flag )  return(0);

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

  for (i=1; i<=mesh->np; i++) {
    ppt = &mesh->point[i];
    ppt->mark = base;
  }
  ++base;
  while ( k != 0 ) {
    pt = &mesh->tria[k];
    if ( !pt->v[0] ) {
      k = pt->nxt;
      continue;
    }
    for (i=0; i<3; i++) {
      ppt = &mesh->point[pt->v[i]];
      if ( ppt->mark != base ) {
        ppt->mark = base;
        s[0] = ppt->c[0];
        s[1] = ppt->c[1];
        s[2] = ppt->c[2];
        if ( ++nbp > MAX_LST-1 ) break;
        listTetraStream(sc,mesh,s,1);
        ppt->flag = 1;
      }
    }
    k = pt->nxt;
  }
  if ( !nbp )  return(0);
  if ( ddebug )  printf("stream start: %d points  ",nbp);
  fprintf(stdout,": %d lines",nbp);
  t = clock() - t;
  fprintf(stdout," %6.2f sec.\n",t/(float)CLOCKS_PER_SEC);

  return(1);
}

int streamRefQuad(pScene sc,pMesh mesh) {
  pMaterial   pm;
  pQuad       pq;
  pPoint      ppt;
  float       s[3];
  int         i,k,nmat,nbp,base;
  time_t      t;

  /* build list */
  base = ++mesh->mark;
  nbp  = 0;
  pq   = &mesh->quad[refitem];

  nmat = matRef(sc,pq->ref);
  /*nmat = !pt->ref ? DEFAULT_MAT : 1+(pt->ref-1)%(sc->par.nbmat-1);*/
  pm = &sc->material[nmat];
  k  = pm->depmat[LQuad];
  if ( !k || pm->flag )  return(0);

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

  for (i=1; i<=mesh->np; i++) {
    ppt = &mesh->point[i];
    ppt->mark = base;
  }
  ++base;
  while ( k != 0 ) {
    pq = &mesh->quad[k];
    if ( !pq->v[0] ) {
      k = pq->nxt;
      continue;
    }
    for (i=0; i<4; i++) {
      ppt = &mesh->point[pq->v[i]];
      if ( ppt->mark != base ) {
        ppt->mark = base;
        s[0] = ppt->c[0];
        s[1] = ppt->c[1];
        s[2] = ppt->c[2];
        if ( ++nbp > MAX_LST-1 ) break;
        listHexaStream(sc,mesh,s,1);
        ppt->flag = 1;
      }
    }
    k = pq->nxt;
  }
  if ( !nbp )  return(0);
  if ( ddebug )  printf("stream start: %d points  ",nbp);
  fprintf(stdout,": %d lines",nbp);
  t = clock() - t;
  fprintf(stdout," %6.2f sec.\n",t/(float)CLOCKS_PER_SEC);

  return(1);
}