src/Plugins/MEDReader/IO/Testing/Cxx/TestMedParallelWrite.cxx

   1 // Copyright (C) 2010-2019  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 /*
  21  * TestMedParallelWrite.cxx
  22  *
  23  *  Created on: 11 mai 2011
  24  *      Author: alejandro
  25  */
  26
  27 #define MED_HAVE_MPI
  28
  29 #include <vtkMed.h>
  30 #define MESGERR 1
  31 #include "med_utils.h"
  32
  33 #include <stdlib.h>
  34 #include <string.h>
  35 #include <assert.h>
  36
  37 int main (int argc, char **argv)
  38 {
  39   const char meshname[MED_NAME_SIZE+1] = "2D unstructured mesh";
  40   const med_int spacedim = 2;
  41   const med_int meshdim = 2;
  42   /*                                         12345678901234561234567890123456 */
  43   const char axisname[2*MED_SNAME_SIZE+1] = "x               y               ";
  44   const char unitname[2*MED_SNAME_SIZE+1] = "cm              cm              ";
  45   med_float coordinates[2222];
  46   const med_int nnodes = 1111;
  47
  48   med_int* quadconnectivity;
  49   const med_int nquad4 = 1000;
  50
  51   med_err _ret=-1;
  52
  53   int mpi_size, mpi_rank;
  54   MPI_Comm comm = MPI_COMM_WORLD;
  55   MPI_Info info = MPI_INFO_NULL;
  56
  57   MPI_Init(&argc, &argv);
  58   MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
  59   MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
  60
  61   med_idt  fid;
  62   char    filename[255]="UsesCase_MEDmesh_parallel.med";
  63   /*     SSCRUTE(_filename); */
  64
  65   if (mpi_rank == 0 ) {
  66     printf("mpi_size = %03d\n", mpi_size);
  67
  68     /* Ouverture du fichier en mode non-parallel */
  69     if ((fid = MEDfileOpen(filename, MED_ACC_CREAT)) < 0){
  70       MED_ERR_(_ret, MED_ERR_OPEN,MED_ERR_FILE, filename);
  71     }
  72
  73     /* write a comment in the file */
  74     if (MEDfileCommentWr(fid,"A 2D unstructured mesh : 15 nodes, 12 cells") < 0) {
  75       MESSAGE("ERROR : write file description ...");
  76     }
  77
  78     /* mesh creation : a 2D unstructured mesh */
  79     if (MEDmeshCr(fid, meshname, spacedim, meshdim, MED_UNSTRUCTURED_MESH,
  80       "A 2D unstructured mesh","",MED_SORT_DTIT,MED_CARTESIAN, axisname, unitname) < 0) {
  81       MESSAGE("ERROR : mesh creation ...");
  82     }
  83
  84     /*
  85      * Building the coordinates of a rectangle of 101 points in the Y-axis,
  86      * and 11 in the X-axis
  87      */
  88     for (int j=0; j<11; j++ )
  89       for (int i=0; i<101; i++ )
  90       {
  91       coordinates[j*202+i*2]   = j+1;
  92       coordinates[j*202+i*2+1] = i+1;
  93       }
  94
  95     /* nodes coordinates in a Cartesian axis in full interlace mode
  96         (X1,Y1, X2,Y2, X3,Y3, ...) with no iteration and computation step
  97      */
  98     if (MEDmeshNodeCoordinateWr(fid, meshname, MED_NO_DT, MED_NO_IT, 0.0,
  99               MED_FULL_INTERLACE, nnodes, coordinates) < 0) {
 100       MESSAGE("ERROR : nodes coordinates ...");
 101     }
 102
 103     if ( MEDfileClose( fid ) < 0) {
 104       MED_ERR_(_ret,MED_ERR_CLOSE,MED_ERR_FILE,"");;
 105     }
 106
 107     MPI_Barrier(comm);
 108   } /* End of process ZERO */
 109   else
 110     {
 111     MPI_Barrier(comm);
 112     }
 113
 114   /* Ouverture du fichier en mode parallel */
 115   if ((fid = MEDparFileOpen(filename, MED_ACC_RDWR ,comm, info)) < 0){
 116     MED_ERR_(_ret, MED_ERR_OPEN,MED_ERR_FILE, filename);
 117   }
 118
 119   med_int     nbofentity = nquad4;
 120   med_int     nbofvaluesperentity = 1;
 121   med_int     nbofconstituentpervalue = 4;
 122   med_int     constituentselect = MED_ALL_CONSTITUENT;
 123   med_switch_mode   switchmode = MED_FULL_INTERLACE;
 124   med_storage_mode    storagemode = MED_COMPACT_STMODE;
 125   const char *const   profilename = MED_NO_PROFILE;
 126
 127   /*
 128    * Calculating block sizes
 129    */
 130
 131   int block_size = (100/mpi_size)*10;
 132   med_size    start  = block_size * mpi_rank + 1;
 133   med_size    stride = block_size;
 134   med_size    count  = 1;
 135   med_size    blocksize = block_size;
 136   med_size    lastblocksize = (100 % mpi_size)*10;
 137   if ((mpi_size == mpi_rank+1) && (lastblocksize != 0))
 138     {
 139     blocksize += lastblocksize;
 140     stride    += lastblocksize;
 141     }
 142   lastblocksize = 0;
 143
 144   printf("%03d: block_size = %03d\n", mpi_rank, block_size);
 145   printf("%03d: start = %03d\n", mpi_rank, start);
 146   printf("%03d: stride = %03d\n", mpi_rank, stride);
 147   printf("%03d: count = %03d\n", mpi_rank, count);
 148   printf("%03d: blocksize = %03d\n", mpi_rank, blocksize);
 149   printf("%03d: lastblocksize = %03d\n", mpi_rank, lastblocksize);
 150   med_filter filter = MED_FILTER_INIT;
 151
 152   if ( MEDfilterBlockOfEntityCr( fid,
 153       nbofentity,
 154       nbofvaluesperentity,
 155       nbofconstituentpervalue,
 156       constituentselect,
 157       switchmode,
 158       storagemode,
 159       profilename,
 160       start,
 161       stride,
 162       count,
 163       blocksize,
 164       lastblocksize,
 165       &filter ) < 0 )
 166     {
 167     MESSAGE("ERROR : filter creation ...");
 168     }
 169
 170   // Attention: there is blocksize and block_size and it does not
 171   // represent the same quantity, in case we are in the last
 172   // block they are different, if not it is the same
 173   quadconnectivity = new med_int[blocksize*4];
 174   int shift = mpi_rank*block_size;
 175   printf("%03d: mpi_rank*block_size = %03d\n", mpi_rank, shift);
 176   printf("%03d: block_size = %03d\n", mpi_rank, block_size);
 177   int base = shift + shift / 101;
 178   int c = 0;
 179   for (int i=0; i<blocksize*4; i+=4 )
 180     {
 181     base++;
 182     if ((base%101) == 0)
 183       base++;
 184
 185     quadconnectivity[i]   = base;
 186     quadconnectivity[i+1] = base+1;
 187     quadconnectivity[i+2] = base+102;
 188     quadconnectivity[i+3] = base+101;
 189     c++;
 190     }
 191   printf("%03d: number of written quads = %03d\n", mpi_rank, c);
 192
 193   if (MEDmeshElementConnectivityAdvancedWr(fid, meshname, MED_NO_DT,
 194            MED_NO_IT, 0.0, MED_CELL, MED_QUAD4,
 195            MED_NODAL, &filter, quadconnectivity) < 0) {
 196     MESSAGE("ERROR : quadrangular cells connectivity ...");
 197   }
 198
 199     if ( MEDfilterClose( &filter ) < 0) {
 200       MESSAGE("ERROR : filter closing ...");
 201     }
 202
 203     if ( MEDfileClose( fid ) < 0) {
 204       MED_ERR_(_ret,MED_ERR_CLOSE,MED_ERR_FILE,""); _ret = -1;
 205     }
 206
 207     /* Barrier before writing family ZERO */
 208     MPI_Barrier(comm);
 209
 210     if (mpi_rank == 0 ) {
 211
 212       if ((fid = MEDfileOpen(filename, MED_ACC_RDWR)) < 0){
 213         MED_ERR_(_ret, MED_ERR_OPEN,MED_ERR_FILE, filename);
 214       }
 215
 216       /* create family 0 : by default, all mesh entities family number is 0 */
 217       if (MEDfamilyCr(fid, meshname,MED_NO_NAME, 0, 0, MED_NO_GROUP) < 0) {
 218         MESSAGE("ERROR : family 0 creation ...");
 219       }
 220
 221       const char familyname_root[MED_NAME_SIZE+1] = "PROCESSOR ";
 222       char familyname[MED_NAME_SIZE+1] = " ";
 223       for (int i=1; i<mpi_size+1; i++)
 224         {
 225         snprintf(familyname, sizeof familyname, "%s%d", familyname_root, i);
 226         if (MEDfamilyCr(fid, meshname,familyname, -i, 0, MED_NO_GROUP) < 0) {
 227           MESSAGE("ERROR : family creation ...");
 228           }
 229         printf("%03d: %s\n", mpi_rank, familyname);
 230         }
 231
 232       med_int familynumbers[nquad4];
 233       int l = 1;
 234       for (int i=0; i<nquad4; i++)
 235         {
 236         if ((i > block_size * l - 1) && (l < mpi_size))
 237           {
 238           l++;
 239           }
 240         familynumbers[i] = -l;
 241         }
 242
 243       if (MEDmeshEntityFamilyNumberWr(fid, meshname, MED_NO_DT, MED_NO_IT,
 244                            MED_CELL, MED_QUAD4, nquad4, familynumbers) < 0) {
 245         MESSAGE("ERROR : nodes family numbers ...");
 246       }
 247
 248       /* Write a Profile */
 249       const char profileName[MED_NAME_SIZE+1] = "QUAD4_PROFILE";
 250       const med_int profilesize = 9;
 251       med_int profilearray[9] = {1, 3, 5, 7, 9, 11, 13, 15, 17};
 252       if (MEDprofileWr(fid, profileName, profilesize, profilearray ) < 0) {
 253         MESSAGE("ERROR : nodes family numbers ...");
 254       }
 255
 256       /* write localization for integration points */
 257       const char localizationName[MED_NAME_SIZE+1] = "QUAD4_INTEGRATION_POINTS_4";
 258       const med_float elementcoordinate[6] = {0.0, 0.0,  1.0, 0.0,  0.0,1.0};
 259       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};
 260       const med_float weight[4] = {1.0/8, 1.0/8, 1.0/8, 1.0/8};
 261       med_int spacedim = 2;
 262       med_int nipoint = 4;
 263       if (MEDlocalizationWr(fid, localizationName, MED_QUAD4, spacedim,
 264           elementcoordinate, MED_FULL_INTERLACE,
 265           nipoint, iPointCoordinate, weight,
 266           MED_NO_INTERPOLATION, MED_NO_MESH_SUPPORT) < 0) {
 267         MESSAGE("ERROR : create family of integration points ...");
 268       }
 269
 270       /* Writing a scalar Field on the Quads right here */
 271       const char fieldname[MED_NAME_SIZE+1] = "TEMPERATURE_FIELD";
 272       const med_int ncomponent = 1;
 273       const char componentname[MED_SNAME_SIZE+1] = "TEMPERATURE";
 274       const char componentunit[MED_SNAME_SIZE+1] = "C";
 275
 276       if (MEDfieldCr(fid, fieldname, MED_FLOAT64,
 277                      ncomponent, componentname, componentunit,"",
 278                      meshname) < 0) {
 279         MESSAGE("ERROR : create field");
 280       }
 281
 282       /* write values at cell (QUADS) centers */
 283       med_float quad4values[nquad4];
 284       for (int i=0; i<nquad4; i++)
 285         quad4values[i] = i%100 + 1;
 286
 287       med_float quad4values4[nquad4 * 4];
 288       long int counter = 0;
 289       for (int i=0; i<nquad4; i++)
 290         {
 291         quad4values[i] = i%100 + 1;
 292         for (int j=0; j<4; j++)
 293           {
 294           quad4values4[counter] = quad4values[i];
 295           counter++;
 296           }
 297         }
 298       if (MEDfieldValueWr(fid, fieldname, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL,
 299                          MED_QUAD4, MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
 300                          nquad4, (unsigned char*) quad4values) < 0) {
 301         MESSAGE("ERROR : write field values on MED_QUAD4");
 302       }
 303
 304       const char fieldname2[MED_NAME_SIZE+1] = "TEMPERATURE_FIELD_PGAUSS";
 305       if (MEDfieldCr(fid, fieldname2, MED_FLOAT64,
 306                      ncomponent, componentname, componentunit,"",
 307                      meshname) < 0) {
 308         MESSAGE("ERROR : create field");
 309       }
 310
 311       if (MEDfieldValueWithProfileWr(
 312              fid, fieldname2, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL, MED_QUAD4,
 313              MED_GLOBAL_PFLMODE, MED_NO_PROFILE, localizationName,
 314            MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
 315            nquad4, (unsigned char*) quad4values4) < 0) {
 316         MESSAGE("ERROR : write field values on MED_QUAD4");
 317       }
 318
 319       const char fieldname3[MED_NAME_SIZE+1] = "TEMPERATURE_FIELD_ELNO";
 320       if (MEDfieldCr(fid, fieldname3, MED_FLOAT64,
 321                      ncomponent, componentname, componentunit,"",
 322                      meshname) < 0) {
 323         MESSAGE("ERROR : create field");
 324       }
 325
 326       if (MEDfieldValueWithProfileWr(
 327              fid, fieldname3, MED_NO_DT, MED_NO_IT, 0.0, MED_CELL, MED_QUAD4,
 328              MED_GLOBAL_PFLMODE, MED_NO_PROFILE, MED_GAUSS_ELNO,
 329            MED_FULL_INTERLACE, MED_ALL_CONSTITUENT,
 330            nquad4, (unsigned char*) quad4values4) < 0) {
 331         MESSAGE("ERROR : write field values on MED_QUAD4");
 332       }
 333
 334       if ( MEDfileClose( fid ) < 0) {
 335         MED_ERR_(_ret,MED_ERR_CLOSE,MED_ERR_FILE,"");;
 336       }
 337
 338       printf("File UsesCase_MEDmesh_parallel.med has been generated.\n");
 339     } /* End of process ZERO */
 340
 341   /* MPI_Finalize must be called AFTER MEDclose which may use MPI calls */
 342   MPI_Finalize();
 343 }