Back to index

salome-paravis  6.5.0
TestMedParallelWrite.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2010-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 /*
00021  * TestMedParallelWrite.cxx
00022  *
00023  *  Created on: 11 mai 2011
00024  *      Author: alejandro
00025  */
00026 
00027 #define MED_HAVE_MPI
00028 
00029 #include <vtkMed.h>
00030 #define MESGERR 1
00031 #include "med_utils.h"
00032 
00033 #include <stdlib.h>
00034 #include <string.h>
00035 #include <assert.h>
00036 
00037 int main (int argc, char **argv)
00038 {
00039   const char meshname[MED_NAME_SIZE+1] = "2D unstructured mesh";
00040   const med_int spacedim = 2;
00041   const med_int meshdim = 2;
00042   /*                                         12345678901234561234567890123456 */
00043   const char axisname[2*MED_SNAME_SIZE+1] = "x               y               ";
00044   const char unitname[2*MED_SNAME_SIZE+1] = "cm              cm              ";
00045   med_float coordinates[2222];
00046   const med_int nnodes = 1111;
00047 
00048   med_int* quadconnectivity;
00049   const med_int nquad4 = 1000;
00050 
00051   med_err _ret=-1;
00052 
00053   int mpi_size, mpi_rank;
00054   MPI_Comm comm = MPI_COMM_WORLD;
00055   MPI_Info info = MPI_INFO_NULL;
00056 
00057   MPI_Init(&argc, &argv);
00058   MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
00059   MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
00060 
00061   med_idt  fid;
00062   char    filename[255]="UsesCase_MEDmesh_parallel.med";
00063   /*     SSCRUTE(_filename); */
00064 
00065   if (mpi_rank == 0 ) {
00066     printf("mpi_size = %03d\n", mpi_size);
00067 
00068     /* Ouverture du fichier en mode non-parallel */
00069     if ((fid = MEDfileOpen(filename, MED_ACC_CREAT)) < 0){
00070       MED_ERR_(_ret, MED_ERR_OPEN,MED_ERR_FILE, filename);
00071     }
00072 
00073     /* write a comment in the file */
00074     if (MEDfileCommentWr(fid,"A 2D unstructured mesh : 15 nodes, 12 cells") < 0) {
00075       MESSAGE("ERROR : write file description ...");
00076     }
00077 
00078     /* mesh creation : a 2D unstructured mesh */
00079     if (MEDmeshCr(fid, meshname, spacedim, meshdim, MED_UNSTRUCTURED_MESH,
00080       "A 2D unstructured mesh","",MED_SORT_DTIT,MED_CARTESIAN, axisname, unitname) < 0) {
00081       MESSAGE("ERROR : mesh creation ...");
00082     }
00083 
00084     /*
00085      * Building the coordinates of a rectangle of 101 points in the Y-axis,
00086      * and 11 in the X-axis
00087      */
00088     for (int j=0; j<11; j++ )
00089       for (int i=0; i<101; i++ )
00090       {
00091       coordinates[j*202+i*2]   = j+1;
00092       coordinates[j*202+i*2+1] = i+1;
00093       }
00094 
00095     /* nodes coordinates in a Cartesian axis in full interlace mode
00096         (X1,Y1, X2,Y2, X3,Y3, ...) with no iteration and computation step
00097      */
00098     if (MEDmeshNodeCoordinateWr(fid, meshname, MED_NO_DT, MED_NO_IT, 0.0,
00099               MED_FULL_INTERLACE, nnodes, coordinates) < 0) {
00100       MESSAGE("ERROR : nodes coordinates ...");
00101     }
00102 
00103     if ( MEDfileClose( fid ) < 0) {
00104       MED_ERR_(_ret,MED_ERR_CLOSE,MED_ERR_FILE,"");;
00105     }
00106 
00107     MPI_Barrier(comm);
00108   } /* End of process ZERO */
00109   else
00110     {
00111     MPI_Barrier(comm);
00112     }
00113 
00114   /* Ouverture du fichier en mode parallel */
00115   if ((fid = MEDparFileOpen(filename, MED_ACC_RDWR ,comm, info)) < 0){
00116     MED_ERR_(_ret, MED_ERR_OPEN,MED_ERR_FILE, filename);
00117   }
00118 
00119   med_int     nbofentity = nquad4;
00120   med_int     nbofvaluesperentity = 1;
00121   med_int     nbofconstituentpervalue = 4;
00122   med_int     constituentselect = MED_ALL_CONSTITUENT;
00123   med_switch_mode   switchmode = MED_FULL_INTERLACE;
00124   med_storage_mode    storagemode = MED_COMPACT_STMODE;
00125   const char *const   profilename = MED_NO_PROFILE;
00126 
00127   /*
00128    * Calculating block sizes
00129    */
00130 
00131   int block_size = (100/mpi_size)*10;
00132   med_size    start  = block_size * mpi_rank + 1;
00133   med_size    stride = block_size;
00134   med_size    count  = 1;
00135   med_size    blocksize = block_size;
00136   med_size    lastblocksize = (100 % mpi_size)*10;
00137   if ((mpi_size == mpi_rank+1) && (lastblocksize != 0))
00138     {
00139     blocksize += lastblocksize;
00140     stride    += lastblocksize;
00141     }
00142   lastblocksize = 0;
00143 
00144   printf("%03d: block_size = %03d\n", mpi_rank, block_size);
00145   printf("%03d: start = %03d\n", mpi_rank, start);
00146   printf("%03d: stride = %03d\n", mpi_rank, stride);
00147   printf("%03d: count = %03d\n", mpi_rank, count);
00148   printf("%03d: blocksize = %03d\n", mpi_rank, blocksize);
00149   printf("%03d: lastblocksize = %03d\n", mpi_rank, lastblocksize);
00150   med_filter filter = MED_FILTER_INIT;
00151 
00152   if ( MEDfilterBlockOfEntityCr( fid,
00153       nbofentity,
00154       nbofvaluesperentity,
00155       nbofconstituentpervalue,
00156       constituentselect,
00157       switchmode,
00158       storagemode,
00159       profilename,
00160       start,
00161       stride,
00162       count,
00163       blocksize,
00164       lastblocksize,
00165       &filter ) < 0 )
00166     {
00167     MESSAGE("ERROR : filter creation ...");
00168     }
00169 
00170   // Attention: there is blocksize and block_size and it does not
00171   // represent the same quantity, in case we are in the last
00172   // block they are different, if not it is the same
00173   quadconnectivity = new med_int[blocksize*4];
00174   int shift = mpi_rank*block_size;
00175   printf("%03d: mpi_rank*block_size = %03d\n", mpi_rank, shift);
00176   printf("%03d: block_size = %03d\n", mpi_rank, block_size);
00177   int base = shift + shift / 101;
00178   int c = 0;
00179   for (int i=0; i<blocksize*4; i+=4 )
00180     {
00181     base++;
00182     if ((base%101) == 0)
00183       base++;
00184 
00185     quadconnectivity[i]   = base;
00186     quadconnectivity[i+1] = base+1;
00187     quadconnectivity[i+2] = base+102;
00188     quadconnectivity[i+3] = base+101;
00189     c++;
00190     }
00191   printf("%03d: number of written quads = %03d\n", mpi_rank, c);
00192 
00193   if (MEDmeshElementConnectivityAdvancedWr(fid, meshname, MED_NO_DT,
00194            MED_NO_IT, 0.0, MED_CELL, MED_QUAD4,
00195            MED_NODAL, &filter, quadconnectivity) < 0) {
00196     MESSAGE("ERROR : quadrangular cells connectivity ...");
00197   }
00198 
00199     if ( MEDfilterClose( &filter ) < 0) {
00200       MESSAGE("ERROR : filter closing ...");
00201     }
00202 
00203     if ( MEDfileClose( fid ) < 0) {
00204       MED_ERR_(_ret,MED_ERR_CLOSE,MED_ERR_FILE,""); _ret = -1;
00205     }
00206 
00207     /* Barrier before writing family ZERO */
00208     MPI_Barrier(comm);
00209 
00210     if (mpi_rank == 0 ) {
00211 
00212       if ((fid = MEDfileOpen(filename, MED_ACC_RDWR)) < 0){
00213         MED_ERR_(_ret, MED_ERR_OPEN,MED_ERR_FILE, filename);
00214       }
00215 
00216       /* create family 0 : by default, all mesh entities family number is 0 */
00217       if (MEDfamilyCr(fid, meshname,MED_NO_NAME, 0, 0, MED_NO_GROUP) < 0) {
00218         MESSAGE("ERROR : family 0 creation ...");
00219       }
00220 
00221       const char familyname_root[MED_NAME_SIZE+1] = "PROCESSOR ";
00222       char familyname[MED_NAME_SIZE+1] = " ";
00223       for (int i=1; i<mpi_size+1; i++)
00224         {
00225         snprintf(familyname, sizeof familyname, "%s%d", familyname_root, i);
00226         if (MEDfamilyCr(fid, meshname,familyname, -i, 0, MED_NO_GROUP) < 0) {
00227           MESSAGE("ERROR : family creation ...");
00228           }
00229         printf("%03d: %s\n", mpi_rank, familyname);
00230         }
00231 
00232       med_int familynumbers[nquad4];
00233       int l = 1;
00234       for (int i=0; i<nquad4; i++)
00235         {
00236         if ((i > block_size * l - 1) && (l < mpi_size))
00237           {
00238           l++;
00239           }
00240         familynumbers[i] = -l;
00241         }
00242 
00243       if (MEDmeshEntityFamilyNumberWr(fid, meshname, MED_NO_DT, MED_NO_IT,
00244                            MED_CELL, MED_QUAD4, nquad4, familynumbers) < 0) {
00245         MESSAGE("ERROR : nodes family numbers ...");
00246       }
00247 
00248       /* Write a Profile */
00249       const char profileName[MED_NAME_SIZE+1] = "QUAD4_PROFILE";
00250       const med_int profilesize = 9;
00251       med_int profilearray[9] = {1, 3, 5, 7, 9, 11, 13, 15, 17};
00252       if (MEDprofileWr(fid, profileName, profilesize, profilearray ) < 0) {
00253         MESSAGE("ERROR : nodes family numbers ...");
00254       }
00255 
00256       /* write localization for integration points */
00257       const char localizationName[MED_NAME_SIZE+1] = "QUAD4_INTEGRATION_POINTS_4";
00258       const med_float elementcoordinate[6] = {0.0, 0.0,  1.0, 0.0,  0.0,1.0};
00259       const med_float iPointCoordinate[8] = {1.0/5, 1.0/5,  3.0/5, 1.0/5,  1.0/5, 3.0/5,  1.0/3, 1.0/3};
00260       const med_float weight[4] = {1.0/8, 1.0/8, 1.0/8, 1.0/8};
00261       med_int spacedim = 2;
00262       med_int nipoint = 4;
00263       if (MEDlocalizationWr(fid, localizationName, MED_QUAD4, spacedim,
00264           elementcoordinate, MED_FULL_INTERLACE,
00265           nipoint, iPointCoordinate, weight,
00266           MED_NO_INTERPOLATION, MED_NO_MESH_SUPPORT) < 0) {
00267         MESSAGE("ERROR : create family of integration points ...");
00268       }
00269 
00270       /* Writing a scalar Field on the Quads right here */
00271       const char fieldname[MED_NAME_SIZE+1] = "TEMPERATURE_FIELD";
00272       const med_int ncomponent = 1;
00273       const char componentname[MED_SNAME_SIZE+1] = "TEMPERATURE";
00274       const char componentunit[MED_SNAME_SIZE+1] = "C";
00275 
00276       if (MEDfieldCr(fid, fieldname, MED_FLOAT64,
00277                      ncomponent, componentname, componentunit,"",
00278                      meshname) < 0) {
00279         MESSAGE("ERROR : create field");
00280       }
00281 
00282       /* write values at cell (QUADS) centers */
00283       med_float quad4values[nquad4];
00284       for (int i=0; i<nquad4; i++)
00285         quad4values[i] = i%100 + 1;
00286 
00287       med_float quad4values4[nquad4 * 4];
00288       long int counter = 0;
00289       for (int i=0; i<nquad4; i++)
00290         {
00291         quad4values[i] = i%100 + 1;
00292         for (int j=0; j<4; j++)
00293           {
00294           quad4values4[counter] = quad4values[i];
00295           counter++;
00296           }
00297         }
00298       if (MEDfieldValueWr(fid, fieldname, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL,
00299                          MED_QUAD4, MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
00300                          nquad4, (unsigned char*) quad4values) < 0) {
00301         MESSAGE("ERROR : write field values on MED_QUAD4");
00302       }
00303 
00304       const char fieldname2[MED_NAME_SIZE+1] = "TEMPERATURE_FIELD_PGAUSS";
00305       if (MEDfieldCr(fid, fieldname2, MED_FLOAT64,
00306                      ncomponent, componentname, componentunit,"",
00307                      meshname) < 0) {
00308         MESSAGE("ERROR : create field");
00309       }
00310 
00311       if (MEDfieldValueWithProfileWr(
00312              fid, fieldname2, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL, MED_QUAD4,
00313              MED_GLOBAL_PFLMODE, MED_NO_PROFILE, localizationName,
00314            MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
00315            nquad4, (unsigned char*) quad4values4) < 0) {
00316         MESSAGE("ERROR : write field values on MED_QUAD4");
00317       }
00318 
00319       const char fieldname3[MED_NAME_SIZE+1] = "TEMPERATURE_FIELD_ELNO";
00320       if (MEDfieldCr(fid, fieldname3, MED_FLOAT64,
00321                      ncomponent, componentname, componentunit,"",
00322                      meshname) < 0) {
00323         MESSAGE("ERROR : create field");
00324       }
00325 
00326       if (MEDfieldValueWithProfileWr(
00327              fid, fieldname3, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL, MED_QUAD4,
00328              MED_GLOBAL_PFLMODE, MED_NO_PROFILE, MED_GAUSS_ELNO,
00329            MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
00330            nquad4, (unsigned char*) quad4values4) < 0) {
00331         MESSAGE("ERROR : write field values on MED_QUAD4");
00332       }
00333 
00334       if ( MEDfileClose( fid ) < 0) {
00335         MED_ERR_(_ret,MED_ERR_CLOSE,MED_ERR_FILE,"");;
00336       }
00337 
00338       printf("File UsesCase_MEDmesh_parallel.med has been generated.\n");
00339     } /* End of process ZERO */
00340 
00341   /* MPI_Finalize must be called AFTER MEDclose which may use MPI calls */
00342   MPI_Finalize();
00343 }