ilists.c

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

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

        /* analyze vertices */
        nbpos = nbneg = nbnul = 0;
        for (l=0; l<4; l++) {
          p0  = &mesh->point[pt->v[l]];
          ps0 = &mesh->sol[pt->v[l]];
          /*if ( ps0->bb < sc->iso.val[0] )  ps0->bb = sc->iso.val[0];*/
          
          if ( ps0->bb > iso )      pos[nbpos++] = l;
          else if ( ps0->bb < iso ) neg[nbneg++] = l;
          else                      nbnul++;
        }
        if ( nbneg == 4 || nbpos == 4 ) {
          k = pt->nxt;
          continue;
        }

        if ( nbneg == 2 && nbpos == 2 ) {
          for (l=0; l<4; l++) {
            k1  = neg[tn[l]];
            k2  = pos[tp[l]];
            p0  = &mesh->point[pt->v[k1]];
            p1  = &mesh->point[pt->v[k2]];
            ps0 = &mesh->sol[pt->v[k1]];
            ps1 = &mesh->sol[pt->v[k2]];
            cc = 0.0f;
            if ( fabs(ps1->bb-ps0->bb) > 0.0f )
              cc = (iso-ps0->bb) / (ps1->bb-ps0->bb);
            cx[l] = p0->c[0]+cc*(p1->c[0]-p0->c[0]);
            cy[l] = p0->c[1]+cc*(p1->c[1]-p0->c[1]);
            cz[l] = p0->c[2]+cc*(p1->c[2]-p0->c[2]);
          }

          /* compute face normal */
          ax = cx[1]-cx[0]; ay = cy[1]-cy[0]; az = cz[1]-cz[0];
          bx = cx[2]-cx[0]; by = cy[2]-cy[0]; bz = cz[2]-cz[0];
          n[0] = ay*bz - az*by;
          n[1] = az*bx - ax*bz;
          n[2] = ax*by - ay*bx;
          d = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
          if ( d > 0.0f ) {
            d = 1.0f / sqrt(d);
            n[0] *= d;  
            n[1] *= d;  
            n[2] *= d;
          }
          glNormal3fv(n);
          glVertex3f(cx[0],cy[0],cz[0]);
          glVertex3f(cx[1],cy[1],cz[1]);
          glVertex3f(cx[2],cy[2],cz[2]);
          
	      glNormal3fv(n);
          glVertex3f(cx[0],cy[0],cz[0]);
          glVertex3f(cx[2],cy[2],cz[2]);
          glVertex3f(cx[3],cy[3],cz[3]);

          if ( ddebug ) {
            fprintf(outv,"%f %f %f 0\n",cx[0],cy[0],cz[0]);
            fprintf(outv,"%f %f %f 0\n",cx[1],cy[1],cz[1]);
            fprintf(outv,"%f %f %f 0\n",cx[2],cy[2],cz[2]);
            fprintf(outv,"%f %f %f 0\n",cx[3],cy[3],cz[3]);

            fprintf(outf,"%d %d %d 0\n",nv+1,nv+2,nv+3);
            fprintf(outf,"%d %d %d 0\n",nv+1,nv+3,nv+4);
          }
          nv+= 4;
          nf+= 2;
        }
        else if ( !nbnul ) {
          for (l=0; l<3; l++) {
            k1 = nbneg == 3 ? neg[l] : pos[l];
            k2 = nbneg == 3 ? pos[0] : neg[0];
            p0 = &mesh->point[pt->v[k1]];
            p1 = &mesh->point[pt->v[k2]];
            ps0 = &mesh->sol[pt->v[k1]];
            ps1 = &mesh->sol[pt->v[k2]];
            cc = 0.0f;
            if ( fabs(ps1->bb-ps0->bb) > 0.0f ) 
              cc = (iso-ps0->bb) / (ps1->bb-ps0->bb);
            cx[l] = p0->c[0]+cc*(p1->c[0]-p0->c[0]);
            cy[l] = p0->c[1]+cc*(p1->c[1]-p0->c[1]);
            cz[l] = p0->c[2]+cc*(p1->c[2]-p0->c[2]);
          }
          /* compute face normal */
          ax = cx[1]-cx[0]; ay = cy[1]-cy[0]; az = cz[1]-cz[0];
          bx = cx[2]-cx[0]; by = cy[2]-cy[0]; bz = cz[2]-cz[0];
          n[0] = ay*bz - az*by;
          n[1] = az*bx - ax*bz;
          n[2] = ax*by - ay*bx;
          d = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
          if ( d > 0.0f ) {
            d = 1.0f / sqrt(d);
            n[0] *= d;
            n[1] *= d;
            n[2] *= d;
          }
          glNormal3fv(n);
          glVertex3f(cx[0],cy[0],cz[0]);
          glVertex3f(cx[1],cy[1],cz[1]);
          glVertex3f(cx[2],cy[2],cz[2]);

          if ( ddebug ) {
            fprintf(outv,"%f %f %f 0\n",cx[0],cy[0],cz[0]);
            fprintf(outv,"%f %f %f 0\n",cx[1],cy[1],cz[1]);
            fprintf(outv,"%f %f %f 0\n",cx[2],cy[2],cz[2]);
            fprintf(outf,"%d %d %d 0\n",nv+1,nv+2,nv+3);
          }
          nv += 3;
          nf += 1;
        }
        k = pt->nxt;
      }
    }
  }
  glEnd();
  glDepthMask(GL_TRUE);
  glDisable(GL_BLEND);

  if ( ddebug ) {
    fclose(outv);
    fclose(outf);
  }
  printf("  Vertices %d   Triangles %d\n",nv,nf);
  glEndList();
  return(dlist);
}


int tetraIsoPOVray(pScene sc,pMesh mesh) {
  FILE      *isofil;
  pTetra     pt;
  pPoint     p0,p1;
  pMaterial  pm;
  pSolution  ps0,ps1;
  double     delta;
  float      iso,cx[4],cy[4],cz[4],cc;
  int        m,k,k1,k2,i,l,pos[4],neg[4],nbpos,nbneg,nbnul;
  char      *ptr,data[128];
  static int tn[4] = {0,0,1,1};
  static int tp[4] = {0,1,1,0};

  /* default */
  if ( !mesh->ntet || !mesh->nbb || mesh->typage == 1 )  return(0);
  if ( ddebug ) printf("create isosurfaces POVray\n");
  if ( egal(sc->iso.val[0],sc->iso.val[MAXISO-1]) )  return(0);
  delta = sc->iso.val[MAXISO-1] - sc->iso.val[0];

  strcpy(data,mesh->name);
  ptr = strstr(data,".mesh");
  if ( ptr )  ptr = '\0';
  strcat(data,".pov"); 
  if ( ddebug )  fprintf(stdout,"  Writing POVRay file %s\n",data);
  isofil = fopen(data,"w");
  if ( !isofil )  return(0);

  for (i=MAXISO-1; i>=0; i--) {
    iso = sc->iso.val[i];

    if ( i == MAXISO-1 )  iso -= 0.001*fabs(iso)/delta;
    else if ( i == 0 )    iso += 0.001*fabs(iso)/delta;

    fprintf(isofil,"\n#declare isosurf%d = mesh {\n",i);

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

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

        /* analyze vertices */
        nbpos = nbneg = nbnul = 0;
        for (l=0; l<4; l++) {
          p0  = &mesh->point[pt->v[l]];
          ps0 = &mesh->sol[pt->v[l]];
          
          if ( ps0->bb > iso )      pos[nbpos++] = l;
          else if ( ps0->bb < iso ) neg[nbneg++] = l;
          else                      nbnul++;
        }
        if ( nbneg == 4 || nbpos == 4 ) {
          k = pt->nxt;
          continue;
        }

        if ( nbneg == 2 && nbpos == 2 ) {
          for (l=0; l<4; l++) {
            k1  = neg[tn[l]];
            k2  = pos[tp[l]];
            p0  = &mesh->point[pt->v[k1]];
            p1  = &mesh->point[pt->v[k2]];
            ps0 = &mesh->sol[pt->v[k1]];
            ps1 = &mesh->sol[pt->v[k2]];
            cc = 0.0f;
            if ( fabs(ps1->bb-ps0->bb) > 0.0f )
              cc = (iso-ps0->bb) / (ps1->bb-ps0->bb);
            cx[l] = p0->c[0]+cc*(p1->c[0]-p0->c[0]);
            cy[l] = p0->c[1]+cc*(p1->c[1]-p0->c[1]);
            cz[l] = p0->c[2]+cc*(p1->c[2]-p0->c[2]);
          }

	  fprintf(isofil,"triangle {\n");
	  fprintf(isofil,"  <%f,%f,%f>,\n",
	  cx[0]+mesh->xtra,cy[0]+mesh->ytra,cz[0]+mesh->ztra);
 	  fprintf(isofil,"  <%f,%f,%f>,\n",
	  cx[1]+mesh->xtra,cy[1]+mesh->ytra,cz[1]+mesh->ztra);
	  fprintf(isofil,"  <%f,%f,%f>\n" ,
	  cx[2]+mesh->xtra,cy[2]+mesh->ytra,cz[2]+mesh->ztra);
          fprintf(isofil,"}\n");

	  fprintf(isofil,"triangle {\n");
	  fprintf(isofil,"  <%f,%f,%f>,\n",
	  cx[0]+mesh->xtra,cy[0]+mesh->ytra,cz[0]+mesh->ztra);
 	  fprintf(isofil,"  <%f,%f,%f>,\n",
	  cx[2]+mesh->xtra,cy[2]+mesh->ytra,cz[2]+mesh->ztra);
	  fprintf(isofil,"  <%f,%f,%f>\n" ,
	  cx[3]+mesh->xtra,cy[3]+mesh->ytra,cz[3]+mesh->ztra);
          fprintf(isofil,"}\n");
        }
        else if ( !nbnul ) {
          for (l=0; l<3; l++) {
            k1 = nbneg == 3 ? neg[l] : pos[l];
            k2 = nbneg == 3 ? pos[0] : neg[0];
            p0 = &mesh->point[pt->v[k1]];
            p1 = &mesh->point[pt->v[k2]];
            ps0 = &mesh->sol[pt->v[k1]];
            ps1 = &mesh->sol[pt->v[k2]];
            cc = 0.0f;
            if ( fabs(ps1->bb-ps0->bb) > 0.0f ) 
              cc = (iso-ps0->bb) / (ps1->bb-ps0->bb);
            cx[l] = p0->c[0]+cc*(p1->c[0]-p0->c[0]);
            cy[l] = p0->c[1]+cc*(p1->c[1]-p0->c[1]);
            cz[l] = p0->c[2]+cc*(p1->c[2]-p0->c[2]);
          }
	  fprintf(isofil,"triangle {\n");
	  fprintf(isofil,"  <%f,%f,%f>,\n",
	  cx[0]+mesh->xtra,cy[0]+mesh->ytra,cz[0]+mesh->ztra);
 	  fprintf(isofil,"  <%f,%f,%f>,\n",
	  cx[1]+mesh->xtra,cy[1]+mesh->ytra,cz[1]+mesh->ztra);
	  fprintf(isofil,"  <%f,%f,%f>\n",
	  cx[2]+mesh->xtra,cy[2]+mesh->ytra,cz[2]+mesh->ztra);
          fprintf(isofil,"}\n");
        }
        k = pt->nxt;
      }
    }
    fprintf(isofil,"}\n");
  }

  fclose(isofil);
  return(1);
}