src/MEDOP/tut/medcoupling/testmed_gendata.py

   1 #!/usr/bin/env python
   2 #  -*- coding: iso-8859-1 -*-
   3
   4 # This script illustrates the basic usage of MEDCoupling and MEDLoader
   5 # to generate test data files for various cases of med operation. It
   6 # illustrates also the usage of numpy to specify the values of the
   7 # fields when defined on a cartesian mesh (grid).
   8 # (gboulant - 11/07/2011)
   9
  10 import MEDCoupling as MC
  11 import MEDLoader as ML
  12
  13 import numpy
  14
  15 #
  16 # ===============================================================
  17 # Helper functions to create meshes
  18 # ===============================================================
  19 #
  20
  21 def createGridMesh(meshName, nbCellsX, nbCellsY):
  22     """
  23     The mesh is created using MEDCoupling. The code below creates a
  24     cartesian mesh as a grid with nbCellsX segments in the X direction
  25     and nbCellsY in the Y direction (nb. cells = nbCellsX * nbCellsY)
  26     """
  27     print "Creating grid mesh of size %sx%s"%(nbCellsX, nbCellsY)
  28     cmesh=MC.MEDCouplingCMesh.New();
  29
  30     # Create X coordinates
  31     nbNodesX = nbCellsX+1
  32     stepX = 0.1
  33     arrX = [float(i * stepX) for i in range(nbNodesX)]
  34     coordsX=MC.DataArrayDouble.New()
  35     coordsX.setValues(arrX,nbNodesX,1)
  36
  37     # Create Y coordinates
  38     nbNodesY = nbCellsY+1
  39     stepY = 0.1
  40     arrY=[float(i * stepY) for i in range(nbNodesY)]
  41     coordsY=MC.DataArrayDouble.New()
  42     coordsY.setValues(arrY,nbNodesY,1)
  43
  44     # Create the grid
  45     cmesh.setCoords(coordsX,coordsY)
  46     cmesh.setName(meshName)
  47
  48     return cmesh
  49
  50 def unstructuredMesh(cartesianMesh):
  51     """
  52     Convert the cartesian mesh in unstructured mesh for the need of
  53     write function of MEDLoader
  54     """
  55     print "Creating unstructured mesh from %s"%(cartesianMesh.getName())
  56     umesh=cartesianMesh.buildUnstructured();
  57     umesh.setName(cartesianMesh.getName())
  58     return umesh
  59
  60 #
  61 # ===============================================================
  62 # Creating a cartesian mesh
  63 # ===============================================================
  64 #
  65 # The size is the number of discrete values in a direction, and then
  66 # corresponds to the number of cells in that direction.
  67 size=80
  68 #size=512
  69
  70
  71 # >>>
  72 # WARNING: remember the problem of tics and spaces. The parameter
  73 # "size" is considered to be a number of cells (intervals). The number
  74 # of nodes in that direction is size+1.
  75 # <<<
  76
  77 nbCellsX = size
  78 nbNodesX = nbCellsX+1
  79
  80 nbCellsY = size # The size could be different than the X size
  81 nbNodesY = nbCellsY+1
  82
  83 meshName = "Grid_%sx%s"%(nbCellsX, nbCellsY)
  84 cmesh = createGridMesh(meshName, nbCellsX, nbCellsY)
  85 umesh = unstructuredMesh(cmesh)
  86 medFileName="gendata.med"
  87 ML.MEDLoader.WriteUMesh(medFileName,umesh,True);
  88
  89 #
  90 # ===============================================================
  91 # Creating a scalar field, working with numpy
  92 # ===============================================================
  93 #
  94
  95 def createField(fieldName,gridMesh,
  96                 numpy2Darray,typeOfField=MC.ON_CELLS,
  97                 iteration=0):
  98     """
  99     The number of values for the fields is deduced from the sizes of
 100     the numpy array. If typeOfField is ON_CELLS, the size is considered
 101     as the number of cells, otherwise it's considered as the number of
 102     nodes. In any case, it must be consistent with the dimensions of
 103     the numpy 2D array.
 104     """
 105     print "Creating field %s with iteration=%s"%(fieldName,iteration)
 106
 107     # The sizes are deduced from the numpy array. Note that if
 108     # typeOfField is ON_CELLS, then the size should correspond to the
 109     # number of cells, while if typeOfField is ON_NODES, then the size
 110     # should correspond to the number of nodes
 111     [sizeX,sizeY] = numpy2Darray.shape
 112
 113     # We first have to reshape the 2D numpy array in a 1D vector that
 114     # concatenate all the rows
 115     data=numpy2Darray.reshape(1,sizeX*sizeY)[0]
 116     # Then, we can create a simple list as required by the MEDCoupling
 117     # DataArrayDouble. Note also the usage of float type because
 118     # MEDCoupling works only with real numbers
 119     listdata=list(data)
 120
 121     # Create the field using the list obtained from the numpy array
 122     field = MC.MEDCouplingFieldDouble.New(typeOfField,MC.ONE_TIME);
 123     field.setName(fieldName);
 124     field.setMesh(gridMesh);
 125     field.setIteration(iteration)
 126     field.setTimeValue(float(iteration))
 127
 128     nbComponents=1 # Only one single component for a scalar field
 129     nbCells=sizeX*sizeY
 130     dataArray=MC.DataArrayDouble.New();
 131     dataArray.setValues(listdata,nbCells,nbComponents)
 132     field.setArray(dataArray);
 133
 134     return field
 135
 136 def writeField(fieldName, numpy2Darray,
 137                typeOfField=MC.ON_CELLS,
 138                iteration=0):
 139
 140     field = createField(fieldName, umesh, numpy2Darray,
 141                         typeOfField, iteration)
 142     createFromScratch=False
 143     ML.MEDLoader.WriteField(medFileName,field,createFromScratch)
 144
 145
 146 def createTestNumpy2DArray(sizeX, sizeY):
 147     """
 148     This illustrates how to create a numpy 2D array for input of the
 149     createField function.
 150     """
 151     rows=[]
 152     for irow in range(sizeY):
 153         row = numpy.arange(start = irow*sizeY,
 154                            stop  = irow*sizeY+sizeX,
 155                            step  = 1,
 156                            dtype='float64')
 157         rows.append(row)
 158
 159     numpy2Darray = numpy.vstack(rows)
 160     return numpy2Darray
 161
 162 def createTestFieldOnCells():
 163     # Test field on cells
 164     numpy2Darray = createTestNumpy2DArray(sizeX=nbCellsX, sizeY=nbCellsY)
 165     writeField("FieldOnCells", numpy2Darray,
 166                typeOfField=MC.ON_CELLS)
 167
 168 def createTestFieldOnNodes():
 169     # Test field on nodes
 170     numpy2Darray = createTestNumpy2DArray(sizeX=nbNodesX, sizeY=nbNodesY)
 171     writeField("FieldOnNodes", numpy2Darray,
 172                typeOfField=MC.ON_NODES)
 173
 174
 175 #
 176 # =================================================
 177 # Creating a time series
 178 # =================================================
 179 #
 180
 181 # -------------------------------------------------
 182 # Simple demo of the principles
 183 # -------------------------------------------------
 184
 185 # In these functions, (x,y) are the indexes of the element in the
 186 # numpy array. Note that theses indexes maps the indexes of the
 187 # cartesian mesh.
 188
 189 # A function can be a simple python function ...
 190 def f1(x,y):
 191     z = 10*x
 192     print "x=%s\ny=%s\nz=%s"%(x,y,z)
 193     return z
 194
 195 # ... but also a more sophisticated callable object, for example to
 196 # defines some parameters
 197 class Function(object):
 198     def __init__(self, sizeX, sizeY, param):
 199         self.sizeX = sizeX
 200         self.sizeY = sizeY
 201         self.param = param
 202
 203     def function(self, x,y):
 204         z = self.param*x
 205         print "x=%s\ny=%s\nz=%s"%(x,y,z)
 206         return z
 207
 208     def __call__(self, x,y):
 209         return self.function(x,y)
 210
 211 fOnNodes=Function(sizeX=nbNodesX, sizeY=nbNodesY, param=10)
 212 fOnCells=Function(sizeX=nbCellsX, sizeY=nbCellsY, param=3)
 213
 214 def createFunctionField_01():
 215     sizeX=nbNodesX
 216     sizeY=nbNodesY
 217     typeOfField=MC.ON_NODES
 218     f=fOnNodes
 219     numpy2Darray = numpy.fromfunction(f,(sizeX,sizeY),dtype='float64')
 220     writeField("FieldOnNodesUsingFunc", numpy2Darray,typeOfField)
 221
 222     f=fOnCells
 223     sizeX=nbCellsX
 224     sizeY=nbCellsY
 225     typeOfField=MC.ON_CELLS
 226     numpy2Darray = numpy.fromfunction(f,(sizeX,sizeY),dtype='float64')
 227     writeField("FieldOnCellsUsingFunc", numpy2Darray,typeOfField)
 228
 229
 230 # -------------------------------------------------
 231 # Using the pyfunctions package to generate data
 232 # -------------------------------------------------
 233
 234 def createNumpy2DArrayWithFunc(sizeX, sizeY, function):
 235     """
 236     @function : a callable than can be used as a function of X.
 237     Typically function should be an instance of Function object
 238     defined in pyfunctions.functions.
 239     """
 240
 241     # X coordinates should range between 0 and 1 to use the normalized
 242     # functions. We have to generate sizeX points:
 243     step=1./sizeX
 244     arrX=[float(i * step) for i in range(sizeX)]
 245
 246     values = function(arrX)
 247
 248     # Then on can create the base row for the numpy 2D array
 249     rowX = numpy.array(values)
 250     # and replicate this row along the Y axis
 251     rows=[]
 252     for irow in range(sizeY):
 253         rows.append(rowX)
 254
 255     numpy2Darray = numpy.vstack(rows)
 256     return numpy2Darray
 257
 258 from pyfunctions.functions import FuncStiffPulse
 259 def createNumpy2DArrayWithFuncStiff(sizeX, sizeY):
 260     f=FuncStiffPulse(xlimit=0.3,stiffness=30,nbPeriods=10)
 261     return createNumpy2DArrayWithFunc(sizeX, sizeY, f)
 262
 263 def createFunctionField_02():
 264     sizeX=nbCellsX
 265     sizeY=nbCellsY
 266     typeOfField=MC.ON_CELLS
 267     numpy2Darray = createNumpy2DArrayWithFuncStiff(sizeX,sizeY)
 268     writeField("FieldOnCellsUsingFunc02", numpy2Darray,typeOfField)
 269
 270     sizeX=nbNodesX
 271     sizeY=nbNodesY
 272     typeOfField=MC.ON_NODES
 273     numpy2Darray = createNumpy2DArrayWithFuncStiff(sizeX,sizeY)
 274     writeField("FieldOnNodesUsingFunc02", numpy2Darray,typeOfField)
 275
 276 #
 277 # =================================================
 278 # Functions to create custom fields for MEDOP tests
 279 # =================================================
 280 #
 281 def createTimeSeries():
 282     """
 283     Create a single med file with a single mesh and a field defined on
 284     several time steps (time series).
 285     """
 286     meshName = "Grid_%sx%s"%(nbCellsX, nbCellsY)
 287     cmesh = createGridMesh(meshName, nbCellsX, nbCellsY)
 288     umesh = unstructuredMesh(cmesh)
 289     medFileName="timeseries.med"
 290     ML.MEDLoader.WriteUMesh(medFileName,umesh,True);
 291
 292     sizeX=nbNodesX
 293     sizeY=nbNodesY
 294     typeOfField=MC.ON_NODES
 295
 296     nbIterations=10
 297     pulseStiffNess = 20
 298     pulseNbPeriods = 10
 299     for iteration in range(nbIterations):
 300         xlimit = float(iteration)/float(nbIterations)
 301         f=FuncStiffPulse(xlimit,stiffness=pulseStiffNess,nbPeriods=pulseNbPeriods)
 302         numpy2Darray = createNumpy2DArrayWithFunc(sizeX,sizeY,f)
 303         field = createField("Pulse",umesh,numpy2Darray,typeOfField,iteration)
 304         ML.MEDLoader.WriteField(medFileName,field,False)
 305
 306 from pyfunctions.functions import FuncStiffExp
 307 def createParametrics():
 308     """
 309     Create 2 med files containing each a mesh (identical) and a field
 310     defined on this mesh in each file.
 311     """
 312     meshName = "Grid_%sx%s_01"%(nbCellsX, nbCellsY)
 313     cmesh = createGridMesh(meshName, nbCellsX, nbCellsY)
 314     umesh = unstructuredMesh(cmesh)
 315
 316     sizeX=nbNodesX
 317     sizeY=nbNodesY
 318     typeOfField=MC.ON_NODES
 319
 320     medFileName="parametric_01.med"
 321     ML.MEDLoader.WriteUMesh(medFileName,umesh,True);
 322     f=FuncStiffExp(xlimit=0.3,stiffness=30)
 323     numpy2Darray = createNumpy2DArrayWithFunc(sizeX,sizeY,f)
 324     fieldName = "StiffExp_01"
 325     field = createField(fieldName,umesh, numpy2Darray,typeOfField)
 326     ML.MEDLoader.WriteField(medFileName,field,False)
 327
 328     medFileName="parametric_02.med"
 329     umesh.setName("Grid_%sx%s_02"%(nbCellsX, nbCellsY))
 330     ML.MEDLoader.WriteUMesh(medFileName,umesh,True);
 331     f=FuncStiffExp(xlimit=0.4,stiffness=30)
 332     numpy2Darray = createNumpy2DArrayWithFunc(sizeX,sizeY,f)
 333     fieldName = "StiffExp_02"
 334     field = createField(fieldName,umesh, numpy2Darray,typeOfField)
 335     ML.MEDLoader.WriteField(medFileName,field,False)
 336
 337 def createParametrics_demo():
 338     """
 339     Create 2 med files containing each a mesh (identical) and a field
 340     defined on this mesh in each file.
 341     """
 342     meshName = "mesh1"
 343     cmesh = createGridMesh(meshName, nbCellsX, nbCellsY)
 344     umesh = unstructuredMesh(cmesh)
 345
 346     sizeX=nbNodesX
 347     sizeY=nbNodesY
 348     typeOfField=MC.ON_NODES
 349
 350     listIteration = [0,1,2,3,4]
 351
 352     medFileName="parametric_01.med"
 353     ML.MEDLoader.WriteUMesh(medFileName,umesh,True);
 354     fieldName = "field1"
 355     for iteration in listIteration:
 356         #f=FuncStiffPulse(xlimit=0.3+0.1*iteration,stiffness=10,nbPeriods=5)
 357         f=FuncStiffExp(xlimit=0.3+0.1*iteration,stiffness=10)
 358         numpy2Darray = createNumpy2DArrayWithFunc(sizeX,sizeY,f)
 359         field = createField(fieldName,umesh, numpy2Darray,typeOfField,iteration)
 360         ML.MEDLoader.WriteField(medFileName,field,False)
 361
 362     medFileName="parametric_02.med"
 363     umesh.setName("mesh2")
 364     ML.MEDLoader.WriteUMesh(medFileName,umesh,True);
 365     fieldName = "field2"
 366     for iteration in listIteration:
 367         #f=FuncStiffPulse(xlimit=0.3+0.1*iteration,stiffness=10,nbPeriods=6)
 368         f=FuncStiffExp(xlimit=0.3+0.1*iteration,stiffness=15)
 369         numpy2Darray = createNumpy2DArrayWithFunc(sizeX,sizeY,f)
 370         field = createField(fieldName,umesh, numpy2Darray,typeOfField,iteration)
 371         ML.MEDLoader.WriteField(medFileName,field,False)
 372
 373
 374
 375 #
 376 # =================================================
 377 # Main runner
 378 # =================================================
 379 #
 380 if __name__ == "__main__":
 381     #createTestFieldOnCells()
 382     #createTestFieldOnNodes()
 383     #createFunctionField_01()
 384     #createFunctionField_02()
 385     #createTimeSeries()
 386     createParametrics_demo()