inout.c

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

#include "medit.h"
#include "libmesh5.h"
#include "extern.h"


int loadMesh(pMesh mesh) {
  pPoint      ppt;
  pEdge       pr;
  pTriangle   pt;
  pQuad       pq;
  pTetra      ptet;
  pHexa       ph;
  double      d,dp1,dp2,dp3,dn[3];
  float      *n,fp1,fp2,fp3;
  int         i,ia,ib,inm,ref,is,k,disc,nn,nt,nq;
  char       *ptr,data[256];

  printf("use loadMesh\n");
  /* default */
  strcpy(data,mesh->name);
  ptr = strstr(data,".mesh");
  if ( !ptr ) {
    strcat(data,".meshb");
    if ( !(inm = GmfOpenMesh(data,GmfRead,&mesh->ver,&mesh->dim)) ) {
      ptr = strstr(data,".mesh");
      *ptr = '\0';
      strcat(data,".mesh");
      if ( !(inm = GmfOpenMesh(data,GmfRead,&mesh->ver,&mesh->dim)) ) {
        fprintf(stderr,"  ** %s  NOT FOUND.\n",data);
        return(-1);
      }
    }
  }
  else if (!(inm = GmfOpenMesh(data,GmfRead,&mesh->ver,&mesh->dim)) ) {
    fprintf(stderr,"  ** %s  NOT FOUND.\n",data);
    return(-1);
  }
  if ( !quiet )  fprintf(stdout,"  Reading %s\n",data);

  /* parse keywords */
  mesh->np   = GmfStatKwd(inm,GmfVertices);
  printf("mesh->np=%i\n",mesh->np);
  mesh->nt   = GmfStatKwd(inm,GmfTriangles);
  printf("mesh->np=%i\n",mesh->nt);
  mesh->nq   = GmfStatKwd(inm,GmfQuadrilaterals);
  mesh->ntet = GmfStatKwd(inm,GmfTetrahedra);
  mesh->nhex = GmfStatKwd(inm,GmfHexahedra);
  mesh->nc   = GmfStatKwd(inm,GmfCorners);
  mesh->nr   = GmfStatKwd(inm,GmfRequiredVertices);
  mesh->na   = GmfStatKwd(inm,GmfEdges);
  mesh->nri  = GmfStatKwd(inm,GmfRidges);
  mesh->nre  = GmfStatKwd(inm,GmfRequiredEdges);
  mesh->nvn  = GmfStatKwd(inm,GmfNormals);
  mesh->ntg  = GmfStatKwd(inm,GmfTangents);
  mesh->ne = mesh->nt + mesh->nq + mesh->ntet + mesh->nhex;

  /* check space dimension */
  if ( mesh->dim < 2 || mesh->dim > 3 ) {
	fprintf(stdout,"  ## Wrong dim\n");
    GmfCloseMesh(inm);
    return(-1);
  }
  /* check if vertices and elements found */
  if ( !mesh->np ) {
    fprintf(stdout,"  ## No vertex\n");
    GmfCloseMesh(inm);
    return(-1);
  }

  /* memory allocation for mesh */
  if ( !zaldy1(mesh) ) {
    GmfCloseMesh(inm);
    return(-1);
  }

  /* read mesh vertices */
  GmfGotoKwd(inm,GmfVertices);
  for (k=1; k<=mesh->np; k++) {
    ppt = &mesh->point[k];
    if ( mesh->ver == GmfFloat ) {
      if ( mesh->dim == 2 ) {
        GmfGetLin(inm,GmfVertices,&fp1,&fp2,&ref);
  	    ppt->c[0] = fp1;
  	    ppt->c[1] = fp2;
	    printf("vertices %f %f %i\n",fp1,fp2,ref);
      }
      else {
        GmfGetLin(inm,GmfVertices,&fp1,&fp2,&fp3,&ref);
  	    ppt->c[0] = fp1;
  	    ppt->c[1] = fp2;
	    ppt->c[2] = fp3;
      }
    }
    else {
      if ( mesh->dim == 2 ) {
        GmfGetLin(inm,GmfVertices,&dp1,&dp2,&ref);
        ppt->c[0] = dp1;
        ppt->c[1] = dp2;
      }
      else {
        GmfGetLin(inm,GmfVertices,&dp1,&dp2,&dp3,&ref);
        ppt->c[0] = dp1;
        ppt->c[1] = dp2;
        ppt->c[2] = dp3;
      }
    }
    ppt->ref = ref & 0x7fff;
    ppt->tag = M_UNUSED;
  }

  /* read mesh triangles */
  disc = 0 ;
  GmfGotoKwd(inm,GmfTriangles);
  for (k=1; k<=mesh->nt; k++) {
    pt = &mesh->tria[k];
    GmfGetLin(inm,GmfTriangles,&pt->v[0],&pt->v[1],&pt->v[2],&ref);
    pt->ref  = ref & 0x7fff;
    for (i=0; i<3; i++) {    
      if ( pt->v[i] < 1 || pt->v[i] > mesh->np ) {
        disc++;
    pt->v[0] = 0;
        break;
      }
      else {
        ppt = &mesh->point[pt->v[i]];
        ppt->tag &= ~M_UNUSED;
      }
    }
  }

  /* mesh quadrilaterals */
  GmfGotoKwd(inm,GmfQuadrilaterals);
  for (k=1; k<=mesh->nq; k++) {
    pq = &mesh->quad[k];
    GmfGetLin(inm,GmfQuadrilaterals,&pq->v[0],&pq->v[1],&pq->v[2],&pq->v[3],&ref);
    pq->ref  = ref & 0x7fff;
    for (i=0; i<4; i++) {    
      if ( pq->v[i] < 1 || pq->v[i] > mesh->np ) {
        disc++;
        pq->v[0] = 0;
        break;
      }
      else {
        ppt = &mesh->point[pq->v[i]];
        ppt->tag &= ~M_UNUSED;
      }
    }
  }

  /* mesh tetrahedra */
  GmfGotoKwd(inm,GmfTetrahedra);
  for (k=1; k<=mesh->ntet; k++) {
    ptet = &mesh->tetra[k];
    GmfGetLin(inm,GmfTetrahedra,&ptet->v[0],&ptet->v[1],&ptet->v[2],&ptet->v[3],&ref);
    ptet->ref  = ref & 0x7fff;
    for (i=0; i<4; i++) {    
      if ( ptet->v[i] < 1 || ptet->v[i] > mesh->np ) {
        disc++;
        ptet->v[0] = 0;
        break;
      }
      else {
        ppt = &mesh->point[ptet->v[i]];
        ppt->tag &= ~M_UNUSED;
      }
    }
  }

  /* mesh hexahedra */
  GmfGotoKwd(inm,GmfHexahedra);
  for (k=1; k<=mesh->nhex; k++) {
    ph = &mesh->hexa[k];
    GmfGetLin(inm,GmfHexahedra,&ph->v[0],&ph->v[1],&ph->v[2],&ph->v[3],
              &ph->v[4],&ph->v[5],&ph->v[6],&ph->v[7],&ref);
    ph->ref  = ref & 0x7fff;
    for (i=0; i<8; i++) {    
      if ( ph->v[i] < 1 || ph->v[i] > mesh->np ) {
        disc++;
        ph->v[0] = 0;
        break;
      }
      else {
        ppt = &mesh->point[ph->v[i]];
        ppt->tag &= ~M_UNUSED;
      }
    }
  }

  /* mesh corners */
  GmfGotoKwd(inm,GmfCorners);
  for (k=1; k<=mesh->nc; k++) {
    GmfGetLin(inm,GmfCorners,&is);
    if ( is < 1 || is > mesh->np )
      disc++;
    else {
      ppt = &mesh->point[is];
      ppt->tag |= M_CORNER;
      ppt->tag &= ~M_UNUSED;
    }
  }

  /* required vertices */
  GmfGotoKwd(inm,GmfRequiredVertices);
  for (k=1; k<=mesh->nr; k++) {
    GmfGetLin(inm,GmfRequiredVertices,&is);
    if ( is < 1 || is > mesh->np )
      disc++;
    else {
      ppt = &mesh->point[is];
      ppt->tag |= M_REQUIRED;
      ppt->tag &= ~M_UNUSED;
    }
  }

  /*mesh edges */
  GmfGotoKwd(inm,GmfEdges);
  for (k=1; k<=mesh->na; k++) {
    GmfGetLin(inm,GmfEdges,&ia,&ib,&ref);
    if ( ia < 1 || ia > mesh->np || ib < 1 || ib > mesh->np )
      disc++;
    else {
      pr = &mesh->edge[k];
      pr->v[0] = ia;
      pr->v[1] = ib;
      pr->ref  = ref & 0x7fff;
      pr->tag  = !pr->ref ? M_NOTAG : M_TAG;
      ppt = &mesh->point[ia];
      ppt->tag &= ~M_UNUSED;
      ppt = &mesh->point[ib];
      ppt->tag &= ~M_UNUSED;
    }
  }

  /* mesh ridges */
  GmfGotoKwd(inm,GmfRidges);
  for (k=1; k<=mesh->nri; k++) {
    GmfGetLin(inm,GmfRidges,&is);
    if ( is < 1 || is > mesh->na )
      disc++;
    else {
      pr = &mesh->edge[is];
      pr->tag |= M_RIDGE;
    }
  }

  /* required edges */
  GmfGotoKwd(inm,GmfRequiredEdges);
  for (k=1; k<=mesh->nre; k++) {
    GmfGetLin(inm,GmfRequiredEdges,&is);
    if ( is < 1 || is > mesh->na )
      disc++;
    else {
      pr = &mesh->edge[is];
      pr->tag |= M_REQUIRED;
    }
  }


  /* mesh normals */
  GmfGotoKwd(inm,GmfNormals);
  for (k=1; k<=mesh->nvn; k++) {
    n = &mesh->extra->n[3*(k-1)+1];
    if ( mesh->ver == GmfFloat )
      GmfGetLin(inm,GmfNormals,&n[0],&n[1],&n[2]);
    else {
      GmfGetLin(inm,GmfNormals,&dn[0],&dn[1],&dn[2]);
      n[0] = dn[0];
      n[1] = dn[1];
      n[2] = dn[2];
    }
    d = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
    if ( d > 0.0 ) {
      d = 1.0 / sqrt(d);
      n[0] *= d;
      n[1] *= d;
      n[2] *= d;
    }
  }

  if ( mesh->nvn ) {
    /* normals at vertices */
    mesh->extra->iv = GmfStatKwd(inm,GmfNormalAtVertices);
    mesh->extra->nv = (int*)M_calloc(mesh->np+1,sizeof(int),"inmesh");
    assert(mesh->extra->nv);
    GmfGotoKwd(inm,GmfNormalAtVertices);
    for (k=1; k<=mesh->extra->iv; k++) {
      GmfGetLin(inm,GmfNormalAtVertices,&nn,&is);
      if ( nn < 1 || nn > mesh->np )
        disc++;
      else
    mesh->extra->nv[nn] = is;
    }

    /* normals at triangle vertices */
    mesh->extra->it = GmfStatKwd(inm,GmfNormalAtTriangleVertices);
    mesh->extra->nt = (int*)M_calloc(3*mesh->nt+1,sizeof(int),"inmesh");
    assert(mesh->extra->nt);
    GmfGotoKwd(inm,GmfNormalAtTriangleVertices);
    for (k=1; k<=mesh->extra->it; k++) {
      GmfGetLin(inm,GmfNormalAtTriangleVertices,&nt,&is,&nn);
      if ( nt < 1 || nt > mesh->nt || is < 1 || is > 3 || nn < 1 || nn > mesh->nvn )
        disc++;
      else
    mesh->extra->nt[3*(nt-1)+is] = nn;
    }

    /*normals at quadrilateral vertices */
    mesh->extra->iq = GmfStatKwd(inm,GmfNormalAtQuadrilateralVertices);
    mesh->extra->nq = (int*)M_calloc(4*mesh->nq+1,sizeof(int),"inmesh");
    assert(mesh->extra->nq);
    GmfGotoKwd(inm,GmfNormalAtQuadrilateralVertices);
    for (k=1; k<=mesh->extra->iq; k++) {
      GmfGetLin(inm,GmfNormalAtQuadrilateralVertices,&nq,&is,&nn);
      if ( nq < 1 || nq > mesh->nq || is < 1 || is > 4 || nn < 1 || nn > mesh->nvn )
        disc++;
      else
    mesh->extra->nq[3*(nq-1)+is] = nn;
    }
  }

  /*mesh tangents */
  GmfGotoKwd(inm,GmfTangents);
  for (k=1; k<=mesh->ntg; k++) {
    n = &mesh->extra->t[3*(k-1)+1];
    if ( mesh->ver == GmfFloat )
      GmfGetLin(inm,GmfTangents,&n[0],&n[1],&n[2]);
    else {
      GmfGetLin(inm,GmfTangents,&dn[0],&dn[1],&dn[2]);
      n[0] = dn[0];
      n[1] = dn[1];
      n[2] = dn[2];
    }
    d = n[0]*n[0] + n[1]*n[1] + n[2]*n[2];
    if ( d > 0.0 ) {
      d = 1.0 / sqrt(d);
      n[0] *= d;
      n[1] *= d;
      n[2] *= d;
    }
  }

  if (mesh->ntg ) {
    /*tangents at vertices */
    mesh->extra->jv = GmfStatKwd(inm,GmfTangentAtVertices);
    mesh->extra->tv = (int*)M_calloc(mesh->np+1,sizeof(int),"inmesh");
    assert(mesh->extra->tv);
    GmfGotoKwd(inm,GmfTangentAtVertices);
    for (k=1; k<=mesh->extra->jv; k++) {
      GmfGetLin(inm,GmfTangentAtVertices,&nn,&is);
      if ( nn < 1 || nn > mesh->np )
        disc++;
      else
    mesh->extra->tv[nn] = is;
    }
    
    /* tangent at edge vertices */
    mesh->extra->je = GmfStatKwd(inm,GmfTangentAtEdgeVertices);
    mesh->extra->te = (int*)M_calloc(2*mesh->na+1,sizeof(int),"inmesh");
    assert(mesh->extra->te);
    GmfGotoKwd(inm,GmfTangentAtEdgeVertices);
    for (k=1; k<=mesh->extra->je; k++) {
      GmfGetLin(inm,GmfTangentAtEdgeVertices,&nt,&is,&nn);
      if ( nt < 1 || nt > mesh->np || is < 1 || is > 2 || nn < 1 || nn > mesh->ntg )
        disc++;
      else
    mesh->extra->te[3*(nt-1)+is] = is;
    }
  }

  GmfCloseMesh(inm);
  if ( disc > 0 ) {
    fprintf(stdout,"  ## %d entities discarded\n",disc);
  }
  return(1);
}


/*mark clipped elements */