vector.c

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

    for (m=0; m<sc->par.nbmat; m++) {
      pm = &sc->material[m];
      k  = pm->depmat[LTria];
      if ( !k || pm->flag )  continue;

      while ( k != 0 ) {
        pt = &mesh->tria[k];
        if ( !pt->v[0] ) {
          k = pt->nxt;
          continue;
        }

        scal  = 0.5*sizeTria(mesh,k);
        epsra = EPST * scal;

        for (i=0; i<3; i++) {
          ppt = &mesh->point[pt->v[i]];
          if ( ppt->mark == mesh->mark )  continue;
          ppt->mark = mesh->mark;

          ps0 = &mesh->sol[pt->v[i]];
          dd  = iso = sqrt(ps0->m[0]*ps0->m[0] + ps0->m[1]*ps0->m[1]);
          if ( dd < epsra )  continue;

          /* color = norm */
          if ( iso < sc->iso.val[0] ) 
            iso = sc->iso.val[0];  
          else if ( iso > sc->iso.val[MAXISO-1] )
            iso = sc->iso.val[MAXISO-1];
          for (ia=0; ia<MAXISO-1; ia++)
            if ( iso < sc->iso.val[ia] )  break;
          kc = (iso-sc->iso.val[ia-1]) / (sc->iso.val[ia] - sc->iso.val[ia-1]);
          hsv[0] = sc->iso.col[ia-1]*(1.0-kc)+sc->iso.col[ia]*kc;
          hsvrgb(hsv,rgb);
          glColor3dv(rgb);

          dd    = 1.0 / dd;
          u[0]  = ps0->m[0] * dd;
          u[1]  = ps0->m[1] * dd;
          cp[0] = ppt->c[0];
          cp[1] = ppt->c[1];
          drawVector2D(cp,u,scal);
        }
        k = pt->nxt; 
      }
    }
    glEnd();
  }
  else {
    glBegin(GL_LINES);
    for (m=0; m<sc->par.nbmat; m++) {
      pm = &sc->material[m];
      k  = pm->depmat[LTria];
      if ( !k || pm->flag )  continue;

      while ( k != 0 ) {
        pt = &mesh->tria[k];
        if ( !pt->v[0] ) {
          k = pt->nxt;
          continue;
        }

        scal  = 0.5*sizeTria(mesh,k);
        epsra = EPST * scal;
        cp[0] = cp[1] = 0.0;
        for (l=0; l<3; l++) {
          ppt = &mesh->point[pt->v[l]];
          cp[0] += ppt->c[0];
          cp[1] += ppt->c[1];
        }
        cp[0] /= 3.0;
        cp[1] /= 3.0;

        /* color = norm */
        ps0 = &mesh->sol[k];
        dd  = iso = sqrt(ps0->m[0]*ps0->m[0] + ps0->m[1]*ps0->m[1]);
        if ( dd < epsra ) {
          k = pt->nxt; 
          continue;
        }
        if ( iso < sc->iso.val[0] ) 
          iso = sc->iso.val[0];  
        else if ( iso > sc->iso.val[MAXISO-1] )
          iso = sc->iso.val[MAXISO-1];
        for (ia=0; ia<MAXISO-1; ia++)
          if ( iso < sc->iso.val[ia] )  break;
        kc = (iso-sc->iso.val[ia-1]) / (sc->iso.val[ia] - sc->iso.val[ia-1]);
        hsv[0] = sc->iso.col[ia-1]*(1.0-kc)+sc->iso.col[ia]*kc;
        hsvrgb(hsv,rgb);
        glColor3dv(rgb);
        
        dd   = 1.0 / dd;
        u[0] = ps0->m[0] * dd;
        u[1] = ps0->m[1] * dd;
        drawVector2D(cp,u,scal);
        k = pt->nxt;
      }
    }
    glEnd();
  }

  glLineWidth(1.0);
  glEndList();
  return(dlist);
}


GLuint listQuad2dVector(pMesh mesh) {
  pMaterial   pm;
  pQuad       pq;
  pPoint      ppt;
  pSolution   ps0;
  pScene      sc;
  double      rgb[3],u[2],epsra,iso,kc,dd,scalemin,scalemax,scal;
  float       cp[2];
  GLuint      dlist = 0;
  int         ia,i,k,l,m;
  static double hsv[3] = { 0.0, 1.0, 0.80 };

  /* default */
  if ( !mesh->nq || !mesh->nbb )  return(0);

  sc    = cv.scene[currentScene()];
  if ( egal(sc->iso.val[0],sc->iso.val[MAXISO-1]) )  return(0);
  if ( ddebug ) printf("create vector list\n");

  /* create display list */
  dlist = glGenLists(1);
  glNewList(dlist,GL_COMPILE);
  if ( glGetError() )  return(0);

  scalemin = sc->dmax * SCALV;
  scalemax = 15*scalemin;
  mesh->mark++;
  glLineWidth(3.0);

  if ( mesh->typage == 2 ) {
    glBegin(GL_LINES);
    for (m=0; m<sc->par.nbmat; m++) {
      pm = &sc->material[m];
      k  = pm->depmat[LQuad];
      if ( !k || pm->flag )  continue;

      while ( k != 0 ) {
        pq = &mesh->quad[k];
        if ( !pq->v[0] ) {
          k = pq->nxt;
          continue;
        }

        scal  = 0.5*sizeQuad(mesh,k);
        epsra = EPST * scal;

        for (i=0; i<4; i++) {
          ppt = &mesh->point[pq->v[i]];
          if ( ppt->mark == mesh->mark )  continue;
          ppt->mark = mesh->mark;

          ps0 = &mesh->sol[pq->v[i]];
          dd  = iso = sqrt(ps0->m[0]*ps0->m[0] + ps0->m[1]*ps0->m[1]);
          if ( dd < epsra )  continue;

          /* color = norm */
          if ( iso < sc->iso.val[0] ) 
            iso = sc->iso.val[0];  
          else if ( iso > sc->iso.val[MAXISO-1] )
            iso = sc->iso.val[MAXISO-1];
          for (ia=0; ia<MAXISO-1; ia++)
            if ( iso < sc->iso.val[ia] )  break;
          kc = (iso-sc->iso.val[ia-1]) / (sc->iso.val[ia] - sc->iso.val[ia-1]);
          hsv[0] = sc->iso.col[ia-1]*(1.0-kc)+sc->iso.col[ia]*kc;
          hsvrgb(hsv,rgb);
          glColor3dv(rgb);

          dd    = 1.0 / dd;
          u[0]  = ps0->m[0] * dd;
          u[1]  = ps0->m[1] * dd;
          cp[0] = ppt->c[0];
          cp[1] = ppt->c[1];
          drawVector2D(cp,u,scal);
        }
        k = pq->nxt; 
      }
    }
    glEnd();
  }
  else {
    glBegin(GL_LINES);
    for (m=0; m<sc->par.nbmat; m++) {
      pm = &sc->material[m];
      k  = pm->depmat[LQuad];
      if ( !k || pm->flag )  continue;

      while ( k != 0 ) {
        pq = &mesh->quad[k];
        if ( !pq->v[0] ) {
          k = pq->nxt;
          continue;
        }

        scal  = 0.5*sizeQuad(mesh,k);
        epsra = EPST * scal;
        cp[0] = cp[1] = 0.0;
        for (l=0; l<4; l++) {
          ppt = &mesh->point[pq->v[l]];
          cp[0] += ppt->c[0];
          cp[1] += ppt->c[1];
        }
        cp[0] *= 0.25;
        cp[1] *= 0.25;

        /* color = norm */
        ps0 = &mesh->sol[k];
        dd  = iso = sqrt(ps0->m[0]*ps0->m[0] + ps0->m[1]*ps0->m[1]);
        if ( dd < epsra ) {
          k = pq->nxt; 
          continue;
        }
        if ( iso < sc->iso.val[0] ) 
          iso = sc->iso.val[0];  
        else if ( iso > sc->iso.val[MAXISO-1] )
          iso = sc->iso.val[MAXISO-1];
        for (ia=0; ia<MAXISO-1; ia++)
          if ( iso < sc->iso.val[ia] )  break;
        kc = (iso-sc->iso.val[ia-1]) / (sc->iso.val[ia] - sc->iso.val[ia-1]);
        hsv[0] = sc->iso.col[ia-1]*(1.0-kc)+sc->iso.col[ia]*kc;
        hsvrgb(hsv,rgb);
        glColor3dv(rgb);
        
        dd   = 1.0 / dd;
        u[0] = ps0->m[0] * dd;
        u[1] = ps0->m[1] * dd;
        drawVector2D(cp,u,scal);
        k = pq->nxt;
      }
    }
    glEnd();
  }

  glLineWidth(1.0);
  glEndList();
  return(dlist);
}


GLuint listTria3dVector(pMesh mesh) {
  pMaterial   pm;
  pTriangle   pt;
  pPoint      ppt;
  pSolution   ps0;
  pScene      sc;
  double      rgb[3],u[3],epsra,iso,kc,dd,scalemin,scalemax,scal;
  float       cp[3];
  GLuint      dlist = 0;
  int         ia,i,k,l,m;
  static double hsv[3] = { 0.0f, 1.0f, 0.80f };

  /* default */
  if ( !mesh->nbb )  return(0);

  sc    = cv.scene[currentScene()];
  /*if ( egal(sc->iso.val[0],sc->iso.val[MAXISO-1]) )  return(0);*/
  if ( ddebug ) printf("create vector list\n");

  /* create display list */
  dlist = glGenLists(1);
  glNewList(dlist,GL_COMPILE);
  if ( glGetError() )  return(0);

  scalemin = sc->dmax * SCALV;
  scalemax = 15*scalemin;
  mesh->mark++;
  glLineWidth(2.0);

  if ( mesh->typage == 2 ) {
    glBegin(GL_LINES);
    for (m=0; m<sc->par.nbmat; m++) {
      pm = &sc->material[m];
      k  = pm->depmat[LTria];
      if ( !k || pm->flag )  continue;

      while ( k != 0 ) {
        pt = &mesh->tria[k];
        if ( !pt->v[0] ) {
          k = pt->nxt;
          continue;
        }

        scal  = sizeTria(mesh,k);
        epsra = EPST * scal;

        for (i=0; i<3; i++) {
          ppt = &mesh->point[pt->v[i]];
          if ( ppt->mark == mesh->mark )  continue;
          ppt->mark = mesh->mark;

          ps0 = &mesh->sol[pt->v[i]];
          dd  = iso = sqrt(ps0->m[0]*ps0->m[0] + ps0->m[1]*ps0->m[1]\
              + ps0->m[2]*ps0->m[2]);
          if ( dd < epsra )  continue;

          /* color = norm */
          if ( iso < sc->iso.val[0] ) 
            iso = sc->iso.val[0];  
          else if ( iso > sc->iso.val[MAXISO-1] )
            iso = sc->iso.val[MAXISO-1];
          for (ia=0; ia<MAXISO-1; ia++)
            if ( iso < sc->iso.val[ia] )  break;
          kc = (iso-sc->iso.val[ia-1]) / (sc->iso.val[ia] - sc->iso.val[ia-1]);
          hsv[0] = sc->iso.col[ia-1]*(1.0-kc)+sc->iso.col[ia]*kc;
          hsvrgb(hsv,rgb);
          glColor3dv(rgb);

          dd    = 1.0 / dd;
          u[0]  = ps0->m[0] * dd;
          u[1]  = ps0->m[1] * dd;
          u[2]  = ps0->m[2] * dd;
          cp[0] = ppt->c[0];
          cp[1] = ppt->c[1];
          cp[2] = ppt->c[2];
          drawVector3D(cp,u,scal);
        }
        k = pt->nxt; 
      }
    }
    glEnd();
  }

  else {
    glBegin(GL_LINES);
    for (m=0; m<sc->par.nbmat; m++) {
      pm = &sc->material[m];
      k  = pm->depmat[LTria];
      if ( !k || pm->flag )  continue;

      while ( k != 0 ) {
        pt = &mesh->tria[k];
        if ( !pt->v[0] ) {
          k = pt->nxt;
          continue;
        }

        scal  = 0.5*sizeTria(mesh,k);
        epsra = EPST * scal;
        cp[0] = cp[1] = 0.0;
        for (l=0; l<3; l++) {
          ppt = &mesh->point[pt->v[l]];
          cp[0] += ppt->c[0];
          cp[1] += ppt->c[1];
          cp[2] += ppt->c[2];
        }
        cp[0] /= 3.0;
        cp[1] /= 3.0;
        cp[2] /= 3.0;

        /* color = norm */
        ps0 = &mesh->sol[k];
        dd  = iso = sqrt(ps0->m[0]*ps0->m[0] + ps0->m[1]*ps0->m[1]\
            + ps0->m[2]*ps0->m[2]);
        if ( dd < epsra ) {
          k = pt->nxt; 
          continue;
        }
        if ( iso < sc->iso.val[0] ) 
          iso = sc->iso.val[0];  
        else if ( iso > sc->iso.val[MAXISO-1] )
          iso = sc->iso.val[MAXISO-1];
        for (ia=0; ia<MAXISO-1; ia++)
          if ( iso < sc->iso.val[ia] )  break;
        kc = (iso-sc->iso.val[ia-1]) / (sc->iso.val[ia] - sc->iso.val[ia-1]);
        hsv[0] = sc->iso.col[ia-1]*(1.0-kc)+sc->iso.col[ia]*kc;
        hsvrgb(hsv,rgb);
        glColor3dv(rgb);
        
        dd   = 1.0 / dd;
        u[0] = ps0->m[0]/dd;
        u[1] = ps0->m[1]/dd;
        u[2] = ps0->m[2]/dd;
        drawVector3D(cp,u,scal);
        k = pt->nxt;
      }
    }
    glEnd();
  }

  glLineWidth(1.0);
  glEndList();
  return(dlist);
}