Concepts et fonctionnalités de base pour créer et lire un maillage

Créer un maillage non structuré

Lire un maillage non structuré

Créer un maillage structuré

Lire un maillage structuré

Créer un maillage avec des familles

Lire un maillage avec des familles

Créer un maillage non structuré

La routine MEDmeshCr / mmhcre permet de créer un maillage dans un fichier MED. Un maillage est identifié par son nom qui est une chaîne de MED_NAME_SIZE caractères (comme tous les noms d'objets du modèle). Lors de la création du maillage, il est nécessaire de préciser outre son nom, sa dimension, la dimension de l'espace de calcul, son type (structuré ou non structuré), l'unité des pas de temps et mode de tri de ces pas de temps (utile pour les maillages évolutifs), ainsi que le type de repère des coordonnées. Dans un premier temps nous ne travaillerons que sur des maillages non évolutifs.

A minima il est nécessaire de définir pour un tel maillage les coordonnées de ses noeuds, ainsi que la connectivité de ses mailles (cf. Connectivités de référence).

L'écriture des coordonnées des noeuds d'un maillage non structuré se fait par l'intermédiaire de la routine MEDmeshNodeCoordinateWr / mmhcow. Cette routine permet l'écriture d'un tableau de coordonnées selon un pas de temps et/ou un numéro d'itération donnés. Comme pour d'autres routines prenant en paramètre des données multi-composantes, on peut choisir le mode de stockage de ces données en mémoire -i.e. MED_FULL_INTERLACE (les composantes d'une entité apparaissent avant celles de l'entité suivante : X1Y1Z1X2Y2Z2...), ou MED_NO_INTERLACE (les composantes ni de toutes les entités apparaissent avant les composantes ni+1 (X1X2....Y1Y2....Z1Z2...).

La routine MEDmeshElementConnectivityWr / mmhcyw permet d'écrire les connectivités d'un type géométrique correspondant à un élément d'un maillage selon un pas de temps et/ou un numéro d'itération donnés. L'écriture s'effectue donc type géométrique par type géométrique pour tous les types d'élements : maille (MED_CELL), face (MED_DESCENDING_FACE) ou arête (MED_DESCENDING_EDGE). Une maille est décrite soit par connectivité nodale (MED_NODAL : liste ordonnée des numéros des noeuds), soit par connectivité descendante (MED_DESCENDING : liste ordonnée des numéros des noeuds, arêtes ou faces selon la dimension de la maille considérée). A noter qu'en connectivité nodale, on ne doit écrire que les mailles, les entités faces et arêtes n'étant utilisées que pour la connectivité descendante. Dans le modèle MED, toutes les entités d'un maillage ont une numérotation globale implicite qui débute à 1. Pour les éléments, la numérotation implicite s'obtient en parcourant l'ensemble des types géométriques selon l'ordre de définition du modèle de données. La règle à appliquer pour définir la taille mémoire "T" nécessaire au stockage de la connectivité d'une entité géométrique est la suivante. Dans le cas de la connectivité nodale, on a "T" = nombre de noeuds. Dans le cas de la connectivité descendante, pour une entité de dimension 3 on a "T" = nombre de faces ; pour une entité de dimension 2, on a "T" = nombre d'arêtes, pour une entité de dimension 1, on a "T" = nombre de noeuds.

Même si on ne définit aucun repérage sur les entités d'un maillage, il est nécessaire de créer la famille de numéro zéro qui ne comporte aucun groupe. Par défaut, toutes les entités du maillage, se rattachent à cette famille. Cette création se fait avec la routine MEDfamilyCr / mfacre.

Le premier exemple suivant permet de créer dans un fichier MED, un maillage non structuré 2D.

/*  This file is part of MED.
 *
 *  COPYRIGHT (C) 1999 - 2015  EDF R&D, CEA/DEN
 *  MED is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MED is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with MED.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 *  How to create an unstructured mesh 
 *
 *  Use case 1 : a 2D unstructured mesh with 15 nodes, 8 triangular cells, 4 triangular cells
 */

#include <med.h>
#define MESGERR 1
#include <med_utils.h>

#include <string.h>

int main (int argc, char **argv) {
  med_idt fid;
  const char meshname[MED_NAME_SIZE+1] = "2D unstructured mesh";
  const med_int spacedim = 2;
  const med_int meshdim = 2;
  /*                                         12345678901234561234567890123456 */
  const char axisname[2*MED_SNAME_SIZE+1] = "x               y               ";
  const char unitname[2*MED_SNAME_SIZE+1] = "cm              cm              ";
  const med_float coordinates[30] = { 2.,1.,  7.,1.,  12.,1.,  17.,1.,  22.,1.,
                                      2.,6.,  7.,6.,  12.,6.,  17.,6.,  22.,6.,
                                      2.,11., 7.,11., 12.,11., 17.,11., 22.,11.};
  const med_int nnodes = 15;
  const med_int triaconnectivity[24] = { 1,7,6,   2,7,1,  3,7,2,   8,7,3,   
                                         13,7,8, 12,7,13, 11,7,12, 6,7,11 };
  const med_int ntria3 = 8;
  const med_int quadconnectivity[16] = {3,4,9,8,    4,5,10,9, 
                                        15,14,9,10, 13,8,9,14};
  const med_int nquad4 = 4;
  med_err ret=-1;
  
  /* open MED file */
  fid = MEDfileOpen("UsesCase_MEDmesh_1.med",MED_ACC_CREAT);
  if (fid < 0) {
    MESSAGE("ERROR : file creation ...");
    goto ERROR;
  }

  /* write a comment in the file */
  if (MEDfileCommentWr(fid,"A 2D unstructured mesh : 15 nodes, 12 cells") < 0) {
    MESSAGE("ERROR : write file description ...");
    goto ERROR;
  }
  
  /* mesh creation : a 2D unstructured mesh */
  if (MEDmeshCr(fid, meshname, spacedim, meshdim, MED_UNSTRUCTURED_MESH, 
                "A 2D unstructured mesh","",MED_SORT_DTIT,MED_CARTESIAN, axisname, unitname) < 0) {
    MESSAGE("ERROR : mesh creation ...");
    goto ERROR;
  }

  /* nodes coordinates in a cartesian axis in full interlace mode 
     (X1,Y1, X2,Y2, X3,Y3, ...) with no iteration and computation step 
  */
  if (MEDmeshNodeCoordinateWr(fid, meshname, MED_NO_DT, MED_NO_IT, 0.0,
                              MED_FULL_INTERLACE, nnodes, coordinates) < 0) {
    MESSAGE("ERROR : nodes coordinates ...");
    goto ERROR;
  }

  /* cells connectiviy is defined in nodal mode with no iteration and computation step */
  if (MEDmeshElementConnectivityWr(fid, meshname, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL, MED_TRIA3,
                                   MED_NODAL, MED_FULL_INTERLACE, ntria3, triaconnectivity) < 0) {
    MESSAGE("ERROR : triangular cells connectivity ...");
    goto ERROR;
  }
  if (MEDmeshElementConnectivityWr(fid, meshname, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL, MED_QUAD4,
                                   MED_NODAL, MED_FULL_INTERLACE, nquad4, quadconnectivity) < 0) {
    MESSAGE("ERROR : quadrangular cells connectivity ...");
    goto ERROR;
  }

  /* create family 0 : by default, all mesh entities family number is 0 */
  if (MEDfamilyCr(fid, meshname,MED_NO_NAME, 0, 0, MED_NO_GROUP) < 0) {
    MESSAGE("ERROR : family 0 creation ...");
    goto ERROR;
  }

  ret = 0;
 ERROR :
  
  /* close MED file */
  if (MEDfileClose(fid)  < 0) {
    MESSAGE("ERROR : close file ...");
    return -1;
  }

  return ret;
}

Lire un maillage non structuré

Dans MED, l'accès aux objets du modèle stockés dans un fichier se fait via deux approches possibles : accès par le nom ou via un itérateur. Pour ce qui concerne les maillages, la routine MEDmeshInfoByName / mmhmin permet de lire les informations relatives à un maillage dont on connaît le nom. Les informations lues correspondent à celles écrites par la routine MEDmeshCr / mmhcre.

Il est ensuite nécessaire pour chaque type d'entité, de lire le nombre d'éléments présents dans le maillage. Cette lecture se faite avec la routine MEDmeshnEntity / mmhnme.

Une fois les nombres d'entités connus (noeuds, mailles, et éventuellement faces et arêtes), il est nécessaire de lire les coordonnées des noeuds avec la routine MEDmeshNodeCoordinateRd / mmhcor et la connectivité de chacun des types géométriques d'éléments présents dans le fichier avec la routine MEDmeshElementConnectivityRd / mmhcyr.

L'exemple suivant permet la lecture du maillage créé avec le premier cas d'utilisation.

/*  This file is part of MED.
 *
 *  COPYRIGHT (C) 1999 - 2015  EDF R&D, CEA/DEN
 *  MED is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MED is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with MED.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 *  Use case 2 : read a 2D unstructured mesh with 15 nodes, 8 triangular cells, 4 triangular cells
 */

#include <med.h>
#define MESGERR 1
#include <med_utils.h>

#include <string.h>

int main (int argc, char **argv) {
  med_idt fid;
  const char meshname[MED_NAME_SIZE+1] = "2D unstructured mesh";
  char meshdescription[MED_COMMENT_SIZE+1];
  med_int meshdim;
  med_int spacedim;
  med_sorting_type sortingtype;
  med_int nstep;
  med_mesh_type meshtype;
  med_axis_type axistype;
  char axisname[2*MED_SNAME_SIZE+1];
  char unitname[2*MED_SNAME_SIZE+1];
  char dtunit[MED_SNAME_SIZE+1];
  med_float *coordinates = NULL;
  med_int nnodes = 0;
  med_int *triaconnectivity = NULL;
  med_int ntria3 = 0;
  med_int *quadconnectivity = NULL;
  med_int nquad4 = 0;
  med_bool coordinatechangement;
  med_bool geotransformation;
  int i;
  int ret = -1;

  /* open MED file with READ ONLY access mode */
  fid = MEDfileOpen("UsesCase_MEDmesh_1.med",MED_ACC_RDONLY);
  if (fid < 0) {
    MESSAGE("ERROR : open file in READ ONLY ACCESS mode ...");
    goto ERROR;
  }

  /* 
   * ... we know that the MED file has only one mesh, 
   * a real code working would check ... 
   */

  /* read mesh informations : mesh dimension, space dimension ... */
  if (MEDmeshInfoByName(fid, meshname, &spacedim, &meshdim, &meshtype, meshdescription, 
                        dtunit, &sortingtype, &nstep, &axistype, axisname, unitname) < 0) {
    MESSAGE("ERROR : mesh info ...");
    goto ERROR;
  }

  /* read how many nodes in the mesh */
  if ((nnodes = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
                               MED_COORDINATE, MED_NO_CMODE,&coordinatechangement,
                               &geotransformation)) < 0) {
    MESSAGE("ERROR : number of nodes ...");
    goto ERROR;
  }
  
  /* 
   * ... we know that we only have MED_TRIA3 and MED_QUAD4 in the mesh, 
   * a real code working would check all MED geometry cell types ... 
   */

  /* read how many triangular cells in the mesh */
  if ((ntria3 = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,MED_TRIA3,
                               MED_CONNECTIVITY, MED_NODAL,&coordinatechangement,
                               &geotransformation)) < 0) {
    MESSAGE("ERROR : number of MED_TRIA3 ...");
    goto ERROR;
  }

  /* read how many quadrangular cells in the mesh */
  if ((nquad4 = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,MED_QUAD4,
                               MED_CONNECTIVITY, MED_NODAL, &coordinatechangement,
                               &geotransformation)) < 0) {
    MESSAGE("ERROR : number of MED_QUAD4 ...");
    goto ERROR;
  }


  /* read mesh nodes coordinates */
  if ((coordinates = (med_float*) malloc(sizeof(med_float)*nnodes*spacedim)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }

  if (MEDmeshNodeCoordinateRd(fid, meshname, MED_NO_DT, MED_NO_IT, MED_FULL_INTERLACE,
                              coordinates) < 0) {
    MESSAGE("ERROR : nodes coordinates ...");
    goto ERROR;
  }

  /* read cells connectivity in the mesh */
  if ((triaconnectivity = (med_int *) malloc(sizeof(med_int)*ntria3*3)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshElementConnectivityRd(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,
                                   MED_TRIA3, MED_NODAL, MED_FULL_INTERLACE, triaconnectivity) < 0) {
    MESSAGE("ERROR : MED_TRIA3 connectivity ...");
    goto ERROR;
  }

  if ((quadconnectivity = (med_int *) malloc(sizeof(med_int)*nquad4*4)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshElementConnectivityRd(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,
                                   MED_QUAD4, MED_NODAL, MED_FULL_INTERLACE, quadconnectivity) < 0) {
    MESSAGE("ERROR : MED_TRIA3 connectivity ...");
    goto ERROR;
  }

  /* 
   * ... we know that the family number of nodes and elements is 0, a real working would check ...
   */
  ret = 0;
 ERROR :

  /* memory desallocation */
  if (coordinates)
    free(coordinates);

  if (triaconnectivity)
    free(triaconnectivity);

  if (quadconnectivity)
    free(quadconnectivity);

  /* close MED file */
  if (MEDfileClose(fid) < 0) {
    MESSAGE("ERROR : close file");             
    ret = -1; 
  } 

  return ret;
}

Une approche plus générique dans la lecture d'un maillage non structuré est donnée dans l'exemple suivant. La routine MEDnMesh / mmhnmh va lire le nombre de maillage dans le fichier. Il s'agit ensuite d'itérer sur ces maillages. La routine MEDmeshInfo / mmhmii permet de lire les informations relatives à un maillage à partir d'un itérateur. La routine MEDmeshnAxis / mmhnax permet de récupérer la dimension de l'espace de calcul en vue de dimensionner la taille des paramètres utiles à la routine MEDmeshInfo.

Comme dans le cas précédent, il est nécessaire de lire les coordonnées des noeuds et la connectivité des mailles du maille (on se place dans le cas d'une connectivité définie en mode nodal). Les routines MEDmeshnEntity / mmhnme et MEDmeshNodeCoordinateRd / mmhcor permettent respectivement de connaître le nombre de noeuds et lire leurs coordonnées. Pour ce qui concerne la lecture de la connectivité des mailles, l'approche proposée est ici plus générique que précédemment. La routine MEDmeshnEntity / mmhnme est d'abord utilisée pour connaître le nombre de type géométriques de mailles présents dans le maillage. Il s'agit ensuite d'itérer sur ce nombre afin de récupérer à chaque itération : le nom et le type géométrique avec la routine MEDmeshEntityInfo / mmheni, le nombre d'élément avec la routine MEDmeshnEntity / mmhnme et le tableau de connectivité avec la routine MEDmeshElementConnectivityRd / mmhcyr.

/*  This file is part of MED.
 *
 *  COPYRIGHT (C) 1999 - 2015  EDF R&D, CEA/DEN
 *  MED is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MED is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with MED.  If not, see <http://www.gnu.org/licenses/>.
 */


/*
 *  Use case 3 : read an unstructured mesh : generic approach
 */


#include <med.h>
#define MESGERR 1
#include <med_utils.h>

#include <string.h>

int main (int argc, char **argv) {
  med_idt fid;
  char meshname[MED_NAME_SIZE+1]="";
  char meshdescription[MED_COMMENT_SIZE+1]="";
  med_int meshdim=0;
  med_int spacedim=0;
  med_sorting_type sortingtype;
  med_int nstep;
  med_mesh_type meshtype;
  med_axis_type axistype;
  char *axisname="";
  char *unitname="";
  char dtunit[MED_SNAME_SIZE+1]="";
  med_float *coordinates = NULL;
  med_int nnodes = 0;
  med_int ngeo = 0;
  med_int nelt=0;
  med_int *connectivity = NULL;
  med_bool coordinatechangement;
  med_bool geotransformation;
  med_int i, j, it, nmesh;
  med_geometry_type geotype;
  med_geometry_type *geotypes = MED_GET_CELL_GEOMETRY_TYPE;
  char geotypename[MED_NAME_SIZE+1];
  int ret=-1;


  /* open MED file with READ ONLY access mode */
  fid = MEDfileOpen("UsesCase_MEDmesh_1.med",MED_ACC_RDONLY);
  if (fid < 0) {
    MESSAGE("ERROR : open file in READ ONLY ACCESS mode ...");
    goto ERROR;
  }

  /* read how many mesh in the file */
  if ((nmesh = MEDnMesh(fid)) < 0) {
    MESSAGE("ERROR : read how many mesh ...");
    goto ERROR;
  }

  for (i=0;i<nmesh;i++) {

    /* read computation space dimension */
    if ((spacedim = MEDmeshnAxis(fid, i+1)) < 0) {
      MESSAGE("ERROR : read computation space dimension ...");
      goto ERROR;
    }

    /* memory allocation */
    if ((axisname  = (char*) malloc(MED_SNAME_SIZE*spacedim+1)) == NULL) {
      MESSAGE("ERROR : memory allocation ...");
      goto ERROR;
    }
    if ((unitname  = (char*) malloc(MED_SNAME_SIZE*spacedim+1)) == NULL) {
      MESSAGE("ERROR : memory allocation ...");
      goto ERROR;
    }
   
    /* read mesh informations : meshname, mesh dimension, mesh type ... */
    if (MEDmeshInfo(fid, i+1, meshname, &spacedim, &meshdim, &meshtype, meshdescription, 
                    dtunit, &sortingtype, &nstep, &axistype, axisname, unitname) < 0) {
      MESSAGE("ERROR : mesh info ...");
      free(axisname); free(unitname);
      goto ERROR;
    }

    free(axisname);
    free(unitname);
    
    /* read how many nodes in the mesh */
    if ((nnodes = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
                                 MED_COORDINATE, MED_NO_CMODE,&coordinatechangement,
                                 &geotransformation)) < 0) {
      MESSAGE("ERROR : number of nodes ...");
    goto ERROR;
    }
    
    /* read mesh nodes coordinates */
    if ((coordinates = (med_float*) malloc(sizeof(med_float)*nnodes*spacedim)) == NULL) {
      MESSAGE("ERROR : memory allocation ...");
      goto ERROR;
    }
    
    if (MEDmeshNodeCoordinateRd(fid, meshname, MED_NO_DT, MED_NO_IT, MED_FULL_INTERLACE,
                                coordinates) < 0) {
      MESSAGE("ERROR : nodes coordinates ...");
      free(coordinates);
      goto ERROR;
    }

    if (coordinates)
      free(coordinates);
    
    /* read number of geometrical types for cells */
    if ((ngeo = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,MED_GEO_ALL,
                               MED_CONNECTIVITY, MED_NODAL,&coordinatechangement,
                               &geotransformation)) < 0) {
      MESSAGE("ERROR : number of cell ...");
      ISCRUTE(geotype);
      goto ERROR;
    }

    for (it=1; it<=ngeo; it++) {

      /* get geometry type */
      if ((nelt = MEDmeshEntityInfo(fid, meshname, MED_NO_DT, MED_NO_IT,MED_CELL,it,
                                   geotypename,&geotype)) < 0) {
        MESSAGE("ERROR : number of cell ...");
        ISCRUTE(it);
        goto ERROR;
      }

      /* how many cells of type geotype ? */
      if ((nelt = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,geotype,
                                 MED_CONNECTIVITY, MED_NODAL,&coordinatechangement,
                                   &geotransformation)) < 0) {
        MESSAGE("ERROR : number of cell ...");
        goto ERROR;
      }
      
      /* read cells connectivity in the mesh */
      if ((connectivity = (med_int *) malloc(sizeof(med_int)*nelt*(geotype%100))) == NULL) {
        MESSAGE("ERROR : memory allocation ...");
        goto ERROR;
      }
      
      if (MEDmeshElementConnectivityRd(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,
                                       geotype, MED_NODAL, MED_FULL_INTERLACE, connectivity) < 0) {
        MESSAGE("ERROR : cellconnectivity ...");
        ISCRUTE(geotype);
        free(connectivity);
        goto ERROR;
      }
      
      /* memory deallocation */
      if (connectivity) {
        free(connectivity);
        connectivity = NULL;
      }


    }
    
  }
  
  /* 
   * ... we know that the family number of nodes and elements is 0, a real working would check ...
   */
  ret=0;
 ERROR:
    
  /* close MED file */
  if (MEDfileClose(fid) < 0) {
    MESSAGE("ERROR : close file");             
    ret=-1; 
  } 

  return ret;
}

Créer un maillage structuré

La création d'un maillage structuré s'effectue avec la routine MEDmeshCr / mmhcre, en précisant qu'il s'agit d'un maillage de type MED_STRUCTURED_MESH. La géométrie d'un maillage structuré est uniquement définie par ses noeuds. La connectivité entre les noeuds est de type (i,j,k) en 3D, (i,j) en 2D et (i) en 1D. On distingue 3 catégories de maillage structuré : les grilles cartésiennes définies dans un repère de coordonnées cartésien (MED_CARTESIAN_GRID), les grilles définies selon un repère de coordonnées polaire (MED_POLAR_GRID), et les grilles curvilignes (MED_CURVILINEAR_GRID).

Dans un maillage structuré, la connaissance des indices d'un noeud en donne la position relative. Nous adopterons le système d'indexation suivant pour ce type de maillage :

indexation_grille.png

Par convention également, le noeud "origine" d'une grille a pour index (1,1,1) en 3D, (1,1) en 2D et (1) en 1D.

Voici par exemple le résultat obtenu pour une grille "5x3".

exemple_grille_cartesienne.png

Dans un maillage structuré, la numérotation des noeuds et des éléments est implicite. Considérons un maillage structuré 2D dont le premier axe comporte "n" points et le second axe "p" points : le noeud d'index (i,j) a pour numéro (i + (j-1) *n) ; l'élément d'index (i,j) (élément porté par le noeud d'index (i,j)) a pour numéro (i + (j-1) * (n-1).

numerotation_grille.png

Si on reporte cette convention de numérotation des noeuds et des éléments sur notre grille "5x3", on obtient :

exemple_numerotation_grille.png



L'écriture des coordonnées d'un maillage structuré correspondant à une grille MED_CARTESIAN_GRID ou MED_POLAR_GRID ne nécessite que l'écriture des indices de coordonnées selon chacun des axes du repère de coordonnées. La routine MEDmeshGridIndexCoordinateWr / mmhgcw permet l'écriture des indices d'un axe donné. Cette routine doit être appelé pour chaque axe du repère de coordonnées. Pour ce qui concerne une grille MED_CURVILINEAR_GRID, il est nécessaire d'écrire toutes les coordonnées des noeuds au même titre que pour un maillage non structuré avec la routine MEDmeshNodeCoordinateWr / mmhcow.

Il est inutile dans un maillage structuré d'écrire la connectivité des mailles dont la définition est implicite. Cependant ces mailles existent bien. Elles sont de type MED_POINT1 pour un maillage 0D, MED_SEG2 pour un maillage 1D, MED_QUAD4 pour un maillage 2D et MED_HEXA8 pour un maillage 3D. Il est possible d'associer à ces mailles tout comme pour les éléments d'un maillage non structuré, les attributs des entités d'un maillage : un numéro (optionnel), un nom (optionnel) de taille MED_SNAME_SIZE caractères, un numéro de famille. La routine MEDmeshEntityNumberWr / mmhenw permet d'écrire un tableau de numéro pour un type d'entité donné. La routine MEDmeshEntityNameWr / mmheaw permet d'écrire un tableau de noms pour un type d'entité donné. L'écriture des noms et numéros des entités d'un maillage est optionnelle. La routine MEDmeshEntityFamilyNumberWr / mmhfnw permet quant à elle l'écriture des numéros de famille. L'écriture des numéros de famille est optionnelle uniquement lorsque toutes les entités du maillage sont rattachés à la famille 0.

/*  This file is part of MED.
 *
 *  COPYRIGHT (C) 1999 - 2015  EDF R&D, CEA/DEN
 *  MED is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MED is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with MED.  If not, see <http://www.gnu.org/licenses/>.
 */


/*
 *  How to create a structured mesh 
 *
 *  Use case 4 :  write a 2D structured mesh (5x3 cartesian grid)
 */

#include <med.h>
#define MESGERR 1
#include <med_utils.h>

#include <string.h>

int main (int argc, char **argv) {
  med_idt fid;
  const char meshname[MED_NAME_SIZE+1] = "2D structured mesh";
  const med_int spacedim = 2;
  const med_int meshdim = 2;
  const char axisname[2*MED_SNAME_SIZE+1] = "x               y               ";
  const char unitname[2*MED_SNAME_SIZE+1] = "cm              cm              ";
  med_int axis, size ;
  const med_float cooXaxis[5] = {1.,2.,3.,4.,5.};
  const med_float cooYaxis[3] = {1.,2.,3.};
  /*                    12345678901234561234567890123456123456789012345612345678901234561234567890123456123456789012345612345678901234561234567890123456 */
  const char cellsnames[8*MED_SNAME_SIZE+1] = "CELL_1          CELL_2          CELL_3          CELL_4          CELL_5          CELL_6          CELL_7          CELL_8          ";
  const med_int nquad4 = 8;
  int ret=-1;

  /* MED file creation */
  fid = MEDfileOpen("UsesCase_MEDmesh_4.med",MED_ACC_CREAT);
  if (fid < 0) {
    MESSAGE("ERROR : file creation ...");
    goto ERROR;
  }

  /* create the structured mesh in the MED file */
  if (MEDmeshCr(fid, meshname, spacedim, meshdim, MED_STRUCTURED_MESH,  
                "A 2D structured mesh","",MED_SORT_DTIT,
                MED_CARTESIAN, axisname, unitname) < 0) {
    MESSAGE("ERROR : mesh creation ...");
    goto ERROR;
  }
  
  /* specify the grid type : MED_CARTESIAN_GRID */
   if (MEDmeshGridTypeWr(fid,meshname, MED_CARTESIAN_GRID) < 0) {
    MESSAGE("ERROR : write grid type ...");
    goto ERROR;
  }
  
   /* write axis "X" and "Y" coordinates */ 
   axis = 1;
   size = 5;
   if (MEDmeshGridIndexCoordinateWr(fid, meshname, MED_NO_DT, MED_NO_IT,0.0, 
                                    axis, size, cooXaxis) < 0) {
     MESSAGE("ERROR : write of axis X coordinates ...");
     goto ERROR;
   }
  axis++;
  size = 3;
  if (MEDmeshGridIndexCoordinateWr(fid, meshname, MED_NO_DT, MED_NO_IT,0.0, 
                                   axis, size, cooYaxis) < 0) {
    MESSAGE("ERROR : write of axis Y coordinates ...");
    goto ERROR;
  }

  /* optionnal : names for nodes or elements */
  /* In this case, a name is given to the cells of the mesh */
  if (MEDmeshEntityNameWr(fid, meshname, MED_NO_DT, MED_NO_IT,
                          MED_CELL, MED_QUAD4, nquad4, 
                          cellsnames) < 0) {
    MESSAGE("ERROR : cells names  ...");
    goto ERROR;
  }
  
  /* create family 0 : by default, all mesh entities family number is 0 */
  if (MEDfamilyCr(fid, meshname,MED_NO_NAME, 0, 0, MED_NO_GROUP) < 0) {
    MESSAGE("ERROR : quadrangular cells connectivity ...");
    goto ERROR;
  }
  
  ret = 0;
 ERROR:

  /* close MED file  */
  if (MEDfileClose(fid) < 0) {
    MESSAGE("ERROR : close file ...");             
    ret = -1; 
  } 

  return ret;
}

Lire un maillage structuré

Une fois le type de maillage identifié avec la routine MEDmeshInfoByName / mmhmin ou MEDmeshInfo / mmhmii, il est nécessaire de lire à quelle catégorie correspond le maillage structuré via la routine MEDmeshGridTypeRd / mmhgtr. Pour les grilles MED_CARTESIAN_GRID ou MED_POLAR_GRID, il est nécessaire de lire les indices de coordonnées selon chaque axes du repère de coordonnées : la routine MEDmeshnEntity / mmhnme permet de lire le nombre d'indices selon un axe donné, la routine MEDmeshGridIndexCoordinateRd / mmhgcr permet la lecture des indices de coordonnées selon les axes de coordonnées du repère. Pour ce qui concerne une grille MED_CURVILINEAR_GRID, il est nécessaire de lire toutes les coordonnées des noeuds au même titre que pour un maillage non structuré avec la routine MEDmeshNodeCoordinateRd / mmhcor.

L'exemple suivant permet la lecture du maillage structuré créé dans le cas d'utilisation précédent.

/*  This file is part of MED.
 *
 *  COPYRIGHT (C) 1999 - 2015  EDF R&D, CEA/DEN
 *  MED is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MED is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with MED.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * Mesh Use case 5 : read a 2D structured mesh
 *                   5x3 cartesian grid
 * 
 */

#include <med.h>
#define MESGERR 1
#include <med_utils.h>

#include <string.h>

int main (int argc, char **argv) {
  med_idt fid;
  const char meshname[MED_NAME_SIZE+1]="2D structured mesh";
  med_int spacedim;
  med_int meshdim;
  char meshdescription[MED_COMMENT_SIZE+1];
  char axisname[2*MED_SNAME_SIZE+1];
  char unitname[2*MED_SNAME_SIZE+1];
  char dtunit[MED_SNAME_SIZE+1];
  med_mesh_type meshtype;
  med_axis_type axistype;
  med_grid_type gridtype;
  med_int axis, size ;
  med_float *cooXaxis = NULL;
  med_float *cooYaxis = NULL;
  med_bool coordinatechangement;
  med_bool geotransformation;
  med_int nstep;
  med_sorting_type sortingtype;
  int j;
  int ret=-1;
  int ncell=0;
  char *cellsname=NULL;

  /* open MED file */
  fid = MEDfileOpen("UsesCase_MEDmesh_4.med",MED_ACC_RDONLY);
  if (fid < 0) {
    MESSAGE("ERROR : open file ...");
    goto ERROR;
  }

  /* read mesh informations : meshname, mesh dimension, space dimension ... */
  if (MEDmeshInfoByName(fid, meshname, &spacedim, &meshdim, &meshtype, meshdescription, 
                        dtunit, &sortingtype, &nstep, &axistype, axisname, unitname) < 0) {
    MESSAGE("ERROR : mesh info ...");
    goto ERROR;
  }
  
  /* read the grid type : MED_CARTESIAN_GRID or MED_CURVILINEAR_GRID */
  if (MEDmeshGridTypeRd(fid, meshname, &gridtype) < 0) {
    MESSAGE("ERROR : read grid type ...");
  }

  /* 
   * ... we know that we the mesh is a cartesian grid, 
   * a real code working would check  ... 
   */

  /* read the axis coordinates (MED_CARTESIAN coordinates system */
  /* X */
  axis = 1;
  if ((size = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, 
                             MED_NODE, MED_NONE, MED_COORDINATE_AXIS1, MED_NO_CMODE,
                             &coordinatechangement, &geotransformation)) < 0) {
    MESSAGE("ERROR : number of coordinates on X axis ...");
    goto ERROR;
  }
  ncell = size-1;

  if ((cooXaxis = (med_float *) malloc(sizeof(med_float)*size)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshGridIndexCoordinateRd(fid, meshname, MED_NO_DT, MED_NO_IT, 
                                   axis, cooXaxis) < 0) {
    MESSAGE("ERROR : read axis X coordinates ...");
    free(cooXaxis);
    goto ERROR;
  }  
  
  free(cooXaxis);

  /* Y */
  axis = 2; 
  if ((size = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, 
                             MED_NODE, MED_NONE, MED_COORDINATE_AXIS2, MED_NO_CMODE,
                             &coordinatechangement, &geotransformation)) < 0) {
    MESSAGE("ERROR : number of coordinates on Y axis ...");
    goto ERROR;
  }
  ncell = ncell * (size-1);

  if ((cooYaxis = (med_float *) malloc(sizeof(med_float)*size)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshGridIndexCoordinateRd(fid, meshname, MED_NO_DT, MED_NO_IT, 
                                   axis, cooYaxis) < 0) {
    MESSAGE("ERROR : read axis Y coordinates ...");
    free(cooYaxis); 
    goto ERROR;
  }

  free(cooYaxis);   

  cellsname = (char *) malloc((sizeof(char))*ncell*MED_SNAME_SIZE+1);
  if (MEDmeshEntityNameRd(fid, meshname, MED_NO_DT, MED_NO_IT, 
                          MED_CELL, MED_QUAD4, cellsname) < 0)  {
    MESSAGE("ERROR : read cells name ...");
    free(cellsname); 
    goto ERROR;
  } 
  free(cellsname);

  ret=0;
 ERROR:
 
  /* close MED file */
  if (MEDfileClose(fid) < 0) {
    MESSAGE("ERROR : close file ...");             
    ret = -1; 
  }

  return ret;
}

Créer un maillage avec des familles

Dans MED, les familles constituent une partition de l'ensemble des entités du maillage : chaque entité du maillage (noeud, maille, face ou arête) appartient à une seule famille. Une famille permet de repérer des zones à particulariser pour le calcul : propriétés physiques, conditions aux limites, extraction des résultats, etc. Un objet famille est identifié par son nom (chaîne de MED_NAME_SIZE caractères) et son numéro. Une famille est décrite par une liste de groupes auxquels les entités de la famille appartiennent. Un groupe est identifié par une chaîne de caractères de MED_LNAME_SIZE caractères.

Une famille donnée ne peut porter que sur un seul type d'entité, une famille est donc soit (exclusif) : une famille de noeuds, une famille d'éléments (mailles/faces/arêtes). Dans un maillage MED, on trouve trois types de familles: la famille de numéro 0 qui ne comporte aucun groupe ; zéro ou plusieurs familles de noeuds dont le numéro doit être strictement positif ; zéro ou plusieurs familles d'éléments (mailles/faces/arêtes) dont le numéro doit être strictement négatif.

La définition de la famille vide de numéro 0 est obligatoire, elle constitue la famille de référence pour tous les noeuds et les éléments qui n'appartiennent à aucun groupe. Une famille de noeuds peut porter le même nom qu'un famille d'éléments. Par contre les familles d'éléments (respectivement de noeuds) doivent toutes avoir des noms différents. La création d'une famille dans un maillage se fait à l'aide de la routine MEDfamilyCr / mfacre.

Le numéro de famille est le lien existant entre la famille considérée et une entité du maillage (noeud, maille, face ou arête). Ce numéro doit être unique dans chaque maillage. La routine MEDmeshEntityFamilyNumberWr / mmhfnw permet l'écriture des numéros de famille pour un type d'entité donné.

Le cas d'utilisation suivant créé dans un maillage non structuré 2D, une famille associée aux noeuds situés sur le bord du maillage. Cette famille porte le numéro 1 et tous les noeuds situés au bord sont rattachés à cette famille. Les numéro de famille des autres noeuds reste à 0.

/*  This file is part of MED.
 *
 *  COPYRIGHT (C) 1999 - 2015  EDF R&D, CEA/DEN
 *  MED is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MED is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with MED.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 *  How to create an unstructured mesh 
 *
 *  Use case 10 : a 2D unstructured mesh with 15 nodes, 8 triangular cells, 4 triangular cells
 *                and families
 */

#include <med.h>
#define MESGERR 1
#include <med_utils.h>

#include <string.h>

int main (int argc, char **argv) {
  med_idt fid;
  const char meshname[MED_NAME_SIZE+1] = "2D unstructured mesh";
  const med_int spacedim = 2;
  const med_int meshdim = 2;
  const char axisname[2*MED_SNAME_SIZE+1] = "x               y               ";
  const char unitname[2*MED_SNAME_SIZE+1] = "cm              cm              ";
  const med_float coordinates[30] = { 2.,1.,  7.,1.,  12.,1.,  17.,1.,  22.,1.,
                                      2.,6.,  7.,6.,  12.,6.,  17.,6.,  22.,6.,
                                      2.,11., 7.,11., 12.,11., 17.,11., 22.,11.};
  const med_int nnodes = 15;
  const med_int triaconnectivity[24] = { 1,7,6,   2,7,1,  3,7,2,   8,7,3,   
                                         13,7,8, 12,7,13, 11,7,12, 6,7,11 };
  const med_int ntria3 = 8;
  const med_int quadconnectivity[16] = {3,4,9,8,    4,5,10,9, 
                                        15,14,9,10, 13,8,9,14};
  const med_int nquad4 = 4;
  const char familyname[MED_NAME_SIZE+1] = "BOUNDARY_VERTICES";
  const char groupname[MED_LNAME_SIZE+1] =  "MESH_BOUNDARY_VERTICES";
  const med_int familynumbers[15] = { 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 1 };
  int ret=-1;
  
  /* open MED file */
  fid = MEDfileOpen("UsesCase_MEDmesh_10.med",MED_ACC_CREAT);
  if (fid < 0) {
    MESSAGE("ERROR : file creation ...");
    goto ERROR;
  }

  /* write a comment in the file */
  if (MEDfileCommentWr(fid,"A 2D unstructured mesh : 15 nodes, 12 cells") < 0) {
    MESSAGE("ERROR : write file description ...");
    goto ERROR;
  }
  
  /* mesh creation : a 2D unstructured mesh */
  if (MEDmeshCr(fid, meshname, spacedim, meshdim, MED_UNSTRUCTURED_MESH, 
                "A 2D unstructured mesh","",MED_SORT_DTIT,MED_CARTESIAN, axisname, unitname) < 0) {
    MESSAGE("ERROR : mesh creation ...");
    goto ERROR;
  }

  /* nodes coordinates in a cartesian axis in full interlace mode 
     (X1,Y1, X2,Y2, X3,Y3, ...) with no iteration and computation step 
  */
  if (MEDmeshNodeCoordinateWr(fid, meshname, MED_NO_DT, MED_NO_IT, 0.0,
                              MED_FULL_INTERLACE, nnodes, coordinates) < 0) {
    MESSAGE("ERROR : nodes coordinates ...");
    goto ERROR;
  }

  /* cells connectiviy is defined in nodal mode */
  if (MEDmeshElementConnectivityWr(fid, meshname, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL, MED_TRIA3,
                                   MED_NODAL, MED_FULL_INTERLACE, ntria3, triaconnectivity) < 0) {
    MESSAGE("ERROR : triangular cells connectivity ...");
    goto ERROR;
  }

  if (MEDmeshElementConnectivityWr(fid, meshname, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL, MED_QUAD4,
                                   MED_NODAL, MED_FULL_INTERLACE, nquad4, quadconnectivity) < 0) {
    MESSAGE("ERROR : quadrangular cells connectivity ...");
    goto ERROR;
  }

  /* create family 0 : by default, all mesh entities family number is 0 */
  if (MEDfamilyCr(fid, meshname,MED_NO_NAME, 0, 0, MED_NO_GROUP) < 0) {
    MESSAGE("ERROR : family creation ...");
    goto ERROR;
  }

  /* create a family for boundary vertices : by convention a nodes family number is > 0, 
     and an element family number is < 0 */
  if (MEDfamilyCr(fid, meshname,familyname, 1, 1, groupname) < 0) {
    MESSAGE("ERROR : family creation ...");
    goto ERROR;
  }

  /* write family number for nodes */
  if (MEDmeshEntityFamilyNumberWr(fid, meshname, MED_NO_DT, MED_NO_IT,
                                  MED_NODE, MED_NONE, nnodes, familynumbers) < 0) {
    MESSAGE("ERROR : nodes family numbers ...");
    goto ERROR;
  }    

  ret=0;
 ERROR:
  
  /* close MED file */
  if (MEDfileClose(fid)  < 0) {
    MESSAGE("ERROR : close file ...");
    ret=-1;
  }

  return ret;
}

Lire un maillage avec des familles

Pour lire un maillage MED de façon générique, il est nécessaire de lire les familles qui outre la famille 0, existent potentiellement dans le maillage. De plus pour chaque type d'entité présent dans le maillage, il est nécessaire de lire le numéro de famille.

La routine MEDnFamily / mfanfa lit le nombre de famille dans un maillage. Il s'agit ensuite d'itérer sur chacune de ces familles. A chaque itération, la routine MEDnFamilyGroup / mfanfg permet de lire le nombre de groupe d'une famille et la routine MEDfamilyInfo / mfafai permet de récupérer les informations relatives à une famille : nom, numéro, liste des groupes. A noter qu'il est également possible d'utiliser la routine MEDfamily23Info / mfaofi pour lire les informations relatives aux familles des fichiers MED 2.3 qui en plus des groupes pouvaient comporter une liste d'attributs entiers (notion devenue obsolète dans MED 3.0), le nombre d'attribut d'une famille pouvant être récupéré avec la routine MEDnFamily23Attribute / mfaona.

Il est possible de lire les numéros de famille pour chaque type d'entité du maillage avec la routine MEDmeshEntityFamilyNumberRd / mmhfnr. Si le tableau de numéro n'est pas présent, cela indique que toutes les entités concernées portent le numéro de famille 0.

Le cas d'utilisation suivant offre un exemple de lecture des familles d'un maillage.

/*  This file is part of MED.
 *
 *  COPYRIGHT (C) 1999 - 2015  EDF R&D, CEA/DEN
 *  MED is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MED is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with MED.  If not, see <http://www.gnu.org/licenses/>.
 */

/*
 *  Use case 11 : read a 2D unstructured mesh with 15 nodes, 8 triangular cells, 4 quadragular cells with
 *  nodes families
 */

#include <med.h>
#define MESGERR 1
#include <med_utils.h>

#include <string.h>

int main (int argc, char **argv) {
  med_idt fid;
  const char meshname[MED_NAME_SIZE+1] = "2D unstructured mesh";
  char meshdescription[MED_COMMENT_SIZE+1]="";
  med_int meshdim;
  med_int spacedim;
  med_sorting_type sortingtype;
  med_int nstep;
  med_mesh_type meshtype;
  med_axis_type axistype;
  char axisname[2*MED_SNAME_SIZE+1]="";
  char unitname[2*MED_SNAME_SIZE+1]="";
  char dtunit[MED_SNAME_SIZE+1];
  med_float *coordinates = NULL;
  med_int nnodes = 0;
  med_int *triaconnectivity = NULL;
  med_int ntria3 = 0;
  med_int *quadconnectivity = NULL;
  med_int nquad4 = 0;
  med_bool coordinatechangement;
  med_bool geotransformation;
  int i;
  med_int nfamily, ngroup;
  med_int familynumber;
  char *groupname=NULL;
  char familyname[MED_NAME_SIZE+1]="";
  med_int *familynumbers = NULL;
  int ret=-1;

  /* open MED file with READ ONLY access mode */
  fid = MEDfileOpen("UsesCase_MEDmesh_10.med",MED_ACC_RDONLY);
  if (fid < 0) {
    MESSAGE("ERROR : open file in READ ONLY ACCESS mode ...");
    goto ERROR;
  }

  /* 
   * ... we know that the MED file has only one mesh, 
   * a real code working would check ... 
   */

  /* read mesh informations : mesh dimension, space dimension ... */
  if (MEDmeshInfoByName(fid, meshname, &spacedim, &meshdim, &meshtype, meshdescription, 
                        dtunit, &sortingtype, &nstep, &axistype, axisname, unitname) < 0) {
    MESSAGE("ERROR : mesh info ...");
    goto ERROR;
  }

  /* read how many nodes in the mesh */
  if ((nnodes = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE,
                               MED_COORDINATE, MED_NO_CMODE,&coordinatechangement,
                               &geotransformation)) < 0) {
    MESSAGE("ERROR : number of nodes ...");
    goto ERROR;
  }
  
  /* 
   * ... we know that we only have MED_TRIA3 and MED_QUAD4 in the mesh, 
   * a real code working would check all MED geometry cell types ... 
   */

  /* read how many triangular cells in the mesh */
  if ((ntria3 = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,MED_TRIA3,
                               MED_CONNECTIVITY, MED_NODAL,&coordinatechangement,
                               &geotransformation)) < 0) {
    MESSAGE("ERROR : number of MED_TRIA3 ...");
    goto ERROR;
  }

  /* read how many quadrangular cells in the mesh */
  if ((nquad4 = MEDmeshnEntity(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,MED_QUAD4,
                               MED_CONNECTIVITY, MED_NODAL, &coordinatechangement,
                               &geotransformation)) < 0) {
    MESSAGE("ERROR : number of MED_QUAD4 ...");
    goto ERROR;
  }

  /* read mesh nodes coordinates */
  if ((coordinates = (med_float*) malloc(sizeof(med_float)*nnodes*spacedim)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }

  if (MEDmeshNodeCoordinateRd(fid, meshname, MED_NO_DT, MED_NO_IT, MED_FULL_INTERLACE,
                              coordinates) < 0) {
    MESSAGE("ERROR : nodes coordinates ...");
    free(coordinates);
    goto ERROR;
  }

  free(coordinates);

  /* read cells connectivity in the mesh */
  if ((triaconnectivity = (med_int *) malloc(sizeof(med_int)*ntria3*3)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshElementConnectivityRd(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,
                                   MED_TRIA3, MED_NODAL, MED_FULL_INTERLACE, triaconnectivity) < 0) {
    MESSAGE("ERROR : MED_TRIA3 connectivity ...");
    free(triaconnectivity);
    goto ERROR;
  }
  free(triaconnectivity);

  if ((quadconnectivity = (med_int *) malloc(sizeof(med_int)*nquad4*4)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshElementConnectivityRd(fid, meshname, MED_NO_DT, MED_NO_IT, MED_CELL,
                                   MED_QUAD4, MED_NODAL, MED_FULL_INTERLACE, quadconnectivity) < 0) {
    MESSAGE("ERROR : MED_TRIA3 connectivity ...");
    free(quadconnectivity);
    goto ERROR;
  }
  free(quadconnectivity);

  /* 
   * read families of entities...
   */
  if ((nfamily = MEDnFamily(fid,meshname)) < 0) {
    MESSAGE("ERROR : read number of family ...");
    goto ERROR;
  }
  
  for (i=0; i<nfamily ; i++) {

    if ((ngroup = MEDnFamilyGroup(fid, meshname, i+1)) < 0) {
      MESSAGE("ERROR : read number of group in a family ...");
      goto ERROR;
    }
    
    if (ngroup > 0) {
      if ((groupname = (char*) malloc(sizeof(char)*MED_LNAME_SIZE*ngroup+1)) == NULL) {
        MESSAGE("ERROR : memory allocation ...");
        goto ERROR;
      }
      
      if (MEDfamilyInfo(fid, meshname, i+1, familyname, &familynumber, groupname) < 0) {
        MESSAGE("ERROR : family info ...");
        free(groupname);
        goto ERROR;
      }
      free(groupname);
    }
    
  }

  /* read family numbers for nodes */
  /* By convention, if there is no numbers in the file, it means that 0 is the family 
     number of all nodes */
  if ((familynumbers = (med_int *) malloc(sizeof(med_int)*nnodes)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshEntityFamilyNumberRd(fid, meshname, MED_NO_DT, MED_NO_IT,
                                  MED_NODE,MED_NONE,familynumbers ) < 0) 
    for (i=0; i<nnodes; i++) *(familynumbers+i) = 0;

  for (i=0; i<nnodes; i++) printf("%d - ", *(familynumbers+i));

  if (familynumbers)
    free (familynumbers);

  /* read family numbers for cells */
  if ((familynumbers = (med_int *) malloc(sizeof(med_int)*ntria3)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshEntityFamilyNumberRd(fid, meshname, MED_NO_DT, MED_NO_IT,
                                  MED_CELL,MED_TRIA3,familynumbers ) < 0) 
    for (i=0; i<ntria3; i++) *(familynumbers+i) = 0;

    free (familynumbers);


  if ((familynumbers = (med_int *) malloc(sizeof(med_int)*nquad4)) == NULL) {
    MESSAGE("ERROR : memory allocation ...");
    goto ERROR;
  }
  if (MEDmeshEntityFamilyNumberRd(fid, meshname, MED_NO_DT, MED_NO_IT,
                                  MED_CELL,MED_QUAD4,familynumbers ) < 0) 
    for (i=0; i<nquad4; i++) *(familynumbers+i) = 0;

  free (familynumbers);

  ret=0;
 ERROR:

  /* close MED file */
  if (MEDfileClose(fid) < 0) {
    MESSAGE("ERROR : close file");             
    ret= -1; 
  } 

  return ret;
}

L'exemple suivant propose une approche similaire avec les routines permettant la lecture des fichiers MED 2.3.

/*  This file is part of MED.
 *
 *  COPYRIGHT (C) 1999 - 2015  EDF R&D, CEA/DEN
 *  MED is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU Lesser General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  MED is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU Lesser General Public License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with MED.  If not, see <http://www.gnu.org/licenses/>.
 */


/******************************************************************************
 * - Nom du fichier : test9.c
 *
 * - Description : lecture des familles d'un maillage MED 
 *
 *****************************************************************************/

#include <med.h>
#define MESGERR 1
#include "med_utils.h"
#include <string.h>

#ifdef DEF_LECT_ECR
#define MODE_ACCES MED_ACC_RDWR
#elif DEF_LECT_AJOUT
#define MODE_ACCES MED_ACC_RDEXT
#else
#define MODE_ACCES MED_ACC_CREAT
#endif

int main (int argc, char **argv)
{
  med_err ret = 0;
  med_idt fid;
  char maa[MED_NAME_SIZE+1];
  med_int mdim;
  med_int nfam;
  int     i,j;
  med_int natt,ngro;
  char *attdes,*gro;
  med_int *attval,*attide;
  char nomfam[MED_NAME_SIZE+1];
  med_int numfam=0, nstep=0, sdim=0;
  char str1[MED_COMMENT_SIZE+1];
  char str2[MED_LNAME_SIZE+1];
  char desc[MED_COMMENT_SIZE+1];
  char dtunit[MED_SNAME_SIZE+1]="";
  char nomcoo[2*MED_SNAME_SIZE+1]="";
  char unicoo[2*MED_SNAME_SIZE+1]="";
  med_mesh_type type;
  med_sorting_type sort;
  med_axis_type rep;

  /* Ouverture du fichier "test8.med" en lecture seule */
  if ((fid = MEDfileOpen("test8.med",MED_ACC_RDONLY)) < 0) {
    MESSAGE("Erreur a l'ouverture du fichier test8.med");
    return -1;
  }

  if ((sdim=MEDmeshnAxis(fid, 1)) <0) {
    MESSAGE("Erreur a la lecture de la dimension de l'espace du maillage :");
    SSCRUTE(maa);
    return -1;
  }

  /* Lecture des infos concernant le premier maillage */
  if ( MEDmeshInfo( fid, 1,  maa, &sdim, &mdim, &type, desc, dtunit, &sort,
                    &nstep,  &rep, nomcoo,unicoo) < 0 ) {
    MESSAGE("Erreur a la lecture des informations sur le maillage : ");SSCRUTE(maa);
    return -1;
  } else {
    printf("Maillage de nom : |%s| , de dimension : "IFORMAT" , et de type %d\n",maa,mdim,type);
    printf("\t -Dimension de l'espace : "IFORMAT"\n",sdim);
    printf("\t -Description du maillage : %s\n",desc);
    printf("\t -Noms des axes : |%s|\n",nomcoo);
    printf("\t -Unités des axes : |%s|\n",unicoo);
    printf("\t -Type de repère : %d\n",rep);
    printf("\t -Nombre d'étape de calcul : "IFORMAT"\n",nstep);
    printf("\t -Unité des dates : |%s|\n",dtunit);
  }


  /* Lecture du nombre de familles */
  if ((nfam = MEDnFamily(fid,maa)) < 0) {
    MESSAGE("Erreur a la lecture du nombre de famille");
    return -1;
  }
  printf("Nombre de familles : "IFORMAT" \n",nfam);

  /* Lecture de chaque famille */
  for (i=0;i<nfam;i++) {

    /* Lecture du nombre de groupe */
    if ((ngro = MEDnFamilyGroup(fid,maa,i+1)) < 0) {
      MESSAGE("Erreur a la lecture du nombre de groupe de la famille d'indice : ");
      ISCRUTE_int(i+1);
      ret = -1;
    }

    /* Lecture du nombre d'attribut */
    if ((natt = MEDnFamily23Attribute(fid,maa,i+1)) < 0) {
      MESSAGE("Erreur a la lecture du nombre d'attribut de la famille d'indice : ");
      ISCRUTE_int(i+1);
      ret = -1;
    }

    if (ret == 0)
      printf("Famille %d a "IFORMAT" attributs et "IFORMAT" groupes \n",i+1,natt,ngro);

    /* Lecture des informations sur la famille */
    if (ret == 0) {
      /* Allocations memoire */
      attide = (med_int*) malloc(sizeof(med_int)*natt);
      attval = (med_int*) malloc(sizeof(med_int)*natt);
      attdes = (char *) malloc(MED_COMMENT_SIZE*natt+1);
      gro = (char*) malloc(MED_LNAME_SIZE*ngro+1);

      if (MEDfamily23Info(fid,maa,i+1,nomfam,attide,attval,attdes,&numfam,gro) < 0) {
        MESSAGE("Erreur a la lecture des informations de la famille d'indice : ");
        ISCRUTE_int(i+1);
        ret = -1;
      }

      if (ret == 0) {
        printf("Famille de nom %s et de numero "IFORMAT" : \n",nomfam,numfam);
        printf("Attributs : \n");
        for (j=0;j<natt;j++) {
          strncpy(str1,attdes+j*MED_COMMENT_SIZE,MED_COMMENT_SIZE);
          str1[MED_COMMENT_SIZE] = '\0';
          printf("ide = "IFORMAT" - val = "IFORMAT" - des = %s\n",*(attide+j),
                 *(attval+j),str1);
        }
        free(attide);
        free(attval);
        free(attdes);

        for (j=0;j<ngro;j++) {
          strncpy(str2,gro+j*MED_LNAME_SIZE,MED_LNAME_SIZE);
          str2[MED_LNAME_SIZE] = '\0';
          printf("gro = %s\n",str2);
              }
        free(gro);
      }
    }
  }

  /* Fermeture du fichier */
  if (MEDfileClose(fid)  < 0) {
    MESSAGE("Erreur a la fermeture du fichier");
    return -1;
  }

  return ret;
}




Généré le Thu Oct 8 14:27:14 2015 pour MED fichier par  doxygen 1.6.1