inout.c

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

static int markPt(pMesh mesh) {
  pTriangle  pt;
  pQuad      pq;
  pPoint     ppt;
  int        i,k,pos,neg,nul;

  for (k=1; k<=mesh->np; k++) {
    ppt = &mesh->point[k];
    if ( ppt->tag & M_UNUSED )  continue;
    ppt->tmp = 0;
  }

  for (k=1; k<=mesh->nt; k++ ) {
    pt = &mesh->tria[k];
    if ( !pt->v[0] )  continue;
    pos = neg = nul = 0;
    for (i=0; i<3; i++) {
      ppt = &mesh->point[pt->v[i]];
      if ( ppt->clip == 2 )       pos++;
      else if ( ppt->clip == 1 )  neg++;
      else                        nul++;
    }
    if ( pos && pos+nul < 4 ) {
      for (i=0; i<3; i++) {
        ppt = &mesh->point[pt->v[i]];
        ppt->tmp = 1;
      }
    }
  }

  for (k=1; k<=mesh->nq; k++ ) {
    pq = &mesh->quad[k];
    if ( !pq->v[0] )  continue;
    pos = neg = nul = 0;
    for (i=0; i<4; i++) {
      ppt = &mesh->point[pq->v[i]];
      if ( ppt->clip == 2 )       pos++;
      else if ( ppt->clip == 1 )  neg++;
      else                        nul++;
    }
    if ( pos && pos+nul < 5 ) {
      for (i=0; i<4; i++) {
        ppt = &mesh->point[pq->v[i]];
        ppt->tmp = 1;
      }
    }
  }

  return(1);
}


/* save (part of) mesh to disk */
int saveMesh(pScene sc,pMesh mesh,char *fileout,ubyte clipon) {
  pPoint     ppt;
  pTriangle  pt;
  pQuad      pq;
  pTetra     ptt;
  pMaterial  pm;
  float      fp1,fp2,fp3;
  int        outm,i,k,m,ver,np,nt,nq,ref;

  ver = GmfFloat;
  if ( !(outm = GmfOpenMesh(fileout,GmfWrite,ver,mesh->dim)) ) {
    fprintf(stderr,"  ** UNABLE TO OPEN %s.\n",fileout);
    return(0);
  }
  fprintf(stdout,"  %%%% %s OPENED\n",fileout);

  /*compact vertices */
  if ( clipon )  
    markPt(mesh);
  np = 0;
  for (k=1; k<=mesh->np; k++) {
    ppt = &mesh->point[k];
    if ( ppt->tag & M_UNUSED )  continue;
    ppt->tmp = ++np;
  }

  GmfSetKwd(outm,GmfVertices,np);
  for (k=1; k<=mesh->np; k++) {
    ppt = &mesh->point[k];
    if ( ppt->tag & M_UNUSED )  continue;
    ref = ppt->ref;
    if ( mesh->dim == 2 ) {
      fp1 = ppt->c[0] + mesh->xtra;
      fp2 = ppt->c[1] + mesh->ytra;
      ref = ppt->ref;
      GmfSetLin(outm,GmfVertices,fp1,fp2,ref);
    }
    else {
      fp1 = ppt->c[0] + mesh->xtra;
      fp2 = ppt->c[1] + mesh->ytra;
      fp3 = ppt->c[2] + mesh->ztra;
      ref = ppt->ref;
      GmfSetLin(outm,GmfVertices,fp1,fp2,fp3,ref);
    }
  }

  /* write triangles */
  nt = 0;
  for (k=1; k<=mesh->nt; k++) {
    pt = &mesh->tria[k];
    if ( !pt->v[0] )  continue;
    m  = matRef(sc,pt->ref);
    pm = &sc->material[m];
    if ( pm->flag )  continue;
    for (i=0; i<3; i++)
      if ( !mesh->point[pt->v[i]].tmp )  break;
    if ( i == 3 )  nt++;
  }
  GmfSetKwd(outm,GmfTriangles,nt);
  for (k=1; k<=mesh->nt; k++) {
    pt = &mesh->tria[k];
    if ( !pt->v[0] )  continue;
    m  = matRef(sc,pt->ref);
    pm = &sc->material[m];
    if ( pm->flag )  continue;
    for (i=0; i<3; i++)
      if ( !mesh->point[pt->v[i]].tmp )  break;
    if ( i < 3 )  continue;
    ref = pt->ref;
    GmfSetLin(outm,GmfTriangles,mesh->point[pt->v[0]].tmp,mesh->point[pt->v[1]].tmp,
              mesh->point[pt->v[2]].tmp,ref);
  }

  /* write quads */
  nq = 0;
  for (k=1; k<=mesh->nq; k++) {
    pq = &mesh->quad[k];
    if ( !pq->v[0] )  continue;
    m  = matRef(sc,pq->ref);
    pm = &sc->material[m];
    if ( pm->flag )  continue;
    for (i=0; i<4; i++)
      if ( !mesh->point[pq->v[i]].tmp )  break;
    if ( i == 4 )  nq++;
  }
  GmfSetKwd(outm,GmfQuadrilaterals,nq);
  for (k=1; k<=mesh->nq; k++) {
    pq = &mesh->quad[k];
    if ( !pq->v[0] )  continue;
    m  = matRef(sc,pq->ref);
    pm = &sc->material[m];
    if ( pm->flag )  continue;
    for (i=0; i<4; i++)
      if ( !mesh->point[pq->v[i]].tmp )  break;
    if ( i < 4 )  continue;
    ref = pq->ref;
    GmfSetLin(outm,GmfQuadrilaterals,mesh->point[pq->v[0]].tmp,mesh->point[pq->v[1]].tmp,
              mesh->point[pq->v[2]].tmp,mesh->point[pq->v[3]].tmp,ref);
  }

  /* write tetrahedra */
  GmfSetKwd(outm,GmfTetrahedra,mesh->ntet);
  for (k=1; k<=mesh->ntet; k++) {
    ptt = &mesh->tetra[k];
    if ( !ptt->v[0] )  continue;
    m  = matRef(sc,ptt->ref);
    pm = &sc->material[m];
    if ( pm->flag )  continue;
    for (i=0; i<4; i++)
      if ( !mesh->point[ptt->v[i]].tmp )  break;
    if ( i < 4 )  continue;
    ref = ptt->ref;
    GmfSetLin(outm,GmfTetrahedra,mesh->point[ptt->v[0]].tmp,mesh->point[ptt->v[1]].tmp,
              mesh->point[ptt->v[2]].tmp,mesh->point[ptt->v[3]].tmp,ref);
  }

  /* write hexahedra */

  if ( !quiet ) {
    fprintf(stdout,"     TOTAL NUMBER OF VERTICES   %8d\n",np);
    fprintf(stdout,"     TOTAL NUMBER OF TRIANGLES  %8d\n",nt);
    fprintf(stdout,"     TOTAL NUMBER OF QUADS      %8d\n",nq);
    fprintf(stdout,"     TOTAL NUMBER OF TETRA      %8d\n",mesh->ntet);
  }

  GmfCloseMesh(outm);
  return(1);
}


/* load solution (metric) */
int loadSol(pMesh mesh,char *filename,int numsol) {
  pSolution    sol;
  double       dbuf[ GmfMaxTyp ];
  float        fbuf[ GmfMaxTyp ];
  double       m[6],lambda[3],eigv[3][3],vp[2][2];
  int          inm,k,i,key,nel,size,type,iord,off,typtab[GmfMaxTyp],ver,dim;
  char        *ptr,data[128];

  strcpy(data,filename);
  ptr = strstr(data,".mesh");
  if ( ptr )  *ptr = '\0';
  strcat(data,".solb");
  if (!(inm = GmfOpenMesh(data,GmfRead,&ver,&dim)) ) {
    ptr  = strstr(data,".sol");
    *ptr = '\0';
    strcat(data,".sol");
    if (!(inm = GmfOpenMesh(data,GmfRead,&ver,&dim)) )
      return(0);
  }
  if ( !quiet )  fprintf(stdout,"  Reading %s\n",data);

  if ( dim != mesh->dim ) {
    fprintf(stderr,"  %%%% Wrong dimension %d.\n",dim);
    GmfCloseMesh(inm);
    return(0);
  }

  nel = GmfStatKwd(inm,GmfSolAtVertices,&type,&size,typtab);
  if ( nel ) {
    if ( nel > mesh->np ) {
      fprintf(stderr,"  %%%% Wrong number: %d Solutions discarded\n",nel-mesh->np);
    }
    mesh->typage = 2;
    key = GmfSolAtVertices;
  }
  else {
    mesh->typage = 1;
    if ( mesh->dim == 2 && mesh->nt ) {
      //if( mesh->nt){
      nel = GmfStatKwd(inm,GmfSolAtTriangles,&type,&size,typtab);
      if ( nel && nel != mesh->nt ) {
        fprintf(stderr,"  %%%% Wrong number %d.\n",nel);
        GmfCloseMesh(inm);
        return(0);
      }
      key = GmfSolAtTriangles;
    }
    else if ( mesh->dim == 2 && mesh->nq ) {
      nel = GmfStatKwd(inm,GmfSolAtQuadrilaterals,&type,&size,typtab);
      if ( nel && nel != mesh->nq ) {
        fprintf(stderr,"  %%%% Wrong number %d.\n",nel);
        GmfCloseMesh(inm);
        return(0);
      }
      key = GmfSolAtQuadrilaterals;
    }
    else if ( mesh->ntet ) {
      nel = GmfStatKwd(inm,GmfSolAtTetrahedra,&type,&size,typtab);
      if ( nel && nel != mesh->ntet ) {
        fprintf(stderr,"  %%%% Wrong number %d.\n",nel);
        GmfCloseMesh(inm);
        return(0);
      }
      key = GmfSolAtTetrahedra;
    }
    else if ( mesh->nhex ) {
      nel = GmfStatKwd(inm,GmfSolAtHexahedra,&type,&size,typtab);
      if ( nel && nel != mesh->nhex ) {
        fprintf(stderr,"  %%%% Wrong number %d.\n",nel);
        GmfCloseMesh(inm);
        return(0);
      }
      key = GmfSolAtHexahedra;
    }
  }
  if ( !nel )  return(1);
  if(ddebug) printf("numsol=%i, type=%i, size=%i\n",numsol,type,size); 
  if ( numsol > type )  numsol = 1;
  numsol--;
  mesh->nbb    = nel;
  mesh->bbmin  =  1.0e20;
  mesh->bbmax  = -1.0e20;
  if ( !zaldy2(mesh) ) {
    GmfCloseMesh(inm);
    return(0);
  }
  sol = mesh->sol;
  sol->dim = dim;
  sol->ver = ver;

  off = 0;
  for (i=0; i<numsol; i++) {
    if(ddebug) printf("typtab[%i]=%i",i,typtab[i]); 
    switch(typtab[i]) {
      case GmfSca:  
        off++; break;
      case GmfVec:
        off += sol->dim; break;
      case GmfSymMat:
        off += sol->dim*(sol->dim+1)/2;  break;
    }
  }
  if(ddebug) printf("typtab[%i]=%i, off%i",i,typtab[i],off); 

  GmfGotoKwd(inm,key);
  if(ddebug) printf("numsol=%i,typtab[i]=%i\n",numsol,typtab[i]); 
  switch(typtab[numsol]) {
    case GmfSca:
      if ( ddebug )  printf("   scalar field\n");
      mesh->nfield = 1;
      for (k=1; k<=nel; k++) {
        if ( sol->ver == GmfFloat )
          GmfGetLin(inm,key,fbuf);
        else {
          GmfGetLin(inm,key,dbuf);
          for (i=0; i<GmfMaxTyp; i++)
            fbuf[i] = dbuf[off+i];
        }
        mesh->sol[k].bb = fbuf[off];
	if(ddebug) printf("valeur donneer %i %f\n",k,fbuf[off]);
	if ( mesh->sol[k].bb < mesh->bbmin )  mesh->bbmin = mesh->sol[k].bb;
        if ( mesh->sol[k].bb > mesh->bbmax )  mesh->bbmax = mesh->sol[k].bb;
      }
      break;

    case GmfVec:
      if ( ddebug )  printf("   vector field\n");
      mesh->nfield = sol->dim;
      for (k=1; k<=nel; k++) {
	mesh->sol[k].bb = 0.0;
        if ( sol->ver == GmfFloat )
          GmfGetLin(inm,key,fbuf);
        else {
          GmfGetLin(inm,key,dbuf);
          for (i=0; i<GmfMaxTyp; i++)
            fbuf[i] = dbuf[off+i];
        }
        for (i=0; i<sol->dim; i++) {
          mesh->sol[k].m[i] = fbuf[off+i];
	  if(ddebug) printf("valeur donner %i composante %i %f\n",k,i,fbuf[off+i]);
          mesh->sol[k].bb  += fbuf[off+i]*fbuf[off+i];
        }
        mesh->sol[k].bb = sqrt(mesh->sol[k].bb);
        if ( mesh->sol[k].bb < mesh->bbmin )  mesh->bbmin = mesh->sol[k].bb;
        if ( mesh->sol[k].bb > mesh->bbmax )  mesh->bbmax = mesh->sol[k].bb;
      }
      break;

    case GmfSymMat:
      if ( ddebug )  printf("   metric field\n");
      mesh->nfield = sol->dim*(sol->dim+1) / 2;
      for (k=1; k<=nel; k++) {
        if ( sol->ver == GmfFloat )
          GmfGetLin(inm,key,fbuf);
        else {
          GmfGetLin(inm,key,dbuf);
          for (i=0; i<GmfMaxTyp; i++)
            fbuf[i] = dbuf[off+i];
        }
        
        if ( sol->dim == 2 ) {
          for (i=0; i<3; i++)  
            mesh->sol[k].m[i] = m[i] = fbuf[off+i];
          iord = eigen2(m,lambda,vp);
          mesh->sol[k].bb = min(lambda[0],lambda[1]);
          if ( mesh->sol[k].bb < mesh->bbmin )  mesh->bbmin = mesh->sol[k].bb;
          if ( mesh->sol[k].bb > mesh->bbmax )  mesh->bbmax = mesh->sol[k].bb;
        }
        else {
          for (i=0; i<6; i++)  
            mesh->sol[k].m[i] = fbuf[off+i];
          mesh->sol[k].m[2] = fbuf[off+3];
          mesh->sol[k].m[3] = fbuf[off+2];
          for (i=0; i<6; i++)  m[i] = mesh->sol[k].m[i];
          iord = eigenv(1,m,lambda,eigv);
          if ( iord ) {
            mesh->sol[k].bb = lambda[0];
            mesh->sol[k].bb = max(mesh->sol[k].bb,lambda[1]);
            mesh->sol[k].bb = max(mesh->sol[k].bb,lambda[2]);
            if ( mesh->sol[k].bb < mesh->bbmin )  mesh->bbmin = mesh->sol[k].bb;
            if ( mesh->sol[k].bb > mesh->bbmax )  mesh->bbmax = mesh->sol[k].bb;
          }
        }
    }
    break;
  }

  GmfCloseMesh(inm);
  return(1);  
}