CDMATH/tests/swig/testFieldCreationAndSave.py

   1 #!/usr/bin/env python3
   2 # -*-coding:utf-8 -*
   3
   4 from math import sqrt
   5 import cdmath
   6 import medcoupling as mc
   7
   8 print("Loading a triangular mesh of a 2D square")
   9 filename = "./meshSquare"
  10 M=cdmath.Mesh(filename+".med", "Mesh_1", 0)
  11
  12 #Extract groups in the mesh
  13 print("Checking boundary group names")
  14 boundaryFaceGroupNames=M.getNameOfFaceGroups()
  15 boundaryNodeGroupNames=M.getNameOfNodeGroups()
  16 print(len(boundaryFaceGroupNames), " Boundary Face Group detected : ", boundaryFaceGroupNames)
  17 print(len(boundaryNodeGroupNames), " Boundary Node Group detected : ", boundaryNodeGroupNames)
  18
  19 assert(len(boundaryFaceGroupNames)==5)
  20 assert(len(boundaryNodeGroupNames)==5)
  21
  22 assert boundaryFaceGroupNames[4]=="Top"
  23 assert boundaryFaceGroupNames[3]=="Right"
  24 assert boundaryFaceGroupNames[2]=="Left"
  25 assert boundaryFaceGroupNames[1]=="Bottom"
  26 assert boundaryFaceGroupNames[0]=="Boundary"
  27
  28 assert boundaryNodeGroupNames[4]=="Top"
  29 assert boundaryNodeGroupNames[3]=="Right"
  30 assert boundaryNodeGroupNames[2]=="Left"
  31 assert boundaryNodeGroupNames[1]=="Bottom"
  32 assert boundaryNodeGroupNames[0]=="Boundary"
  33
  34 #Extract domain sizes
  35 xmin = M.getXMin()
  36 xmax = M.getXMax()
  37 ymin = M.getYMin()
  38 ymax = M.getYMax()
  39
  40 radius = min(xmax-xmin,ymax-ymin)/4
  41 xcentre = (xmax+xmin)/2
  42 ycentre = (ymax+ymin)/2
  43
  44 nbCells = M.getNumberOfCells()
  45 nbNodes = M.getNumberOfNodes()
  46
  47 # Create solid temperature fields
  48 temperature_field_cells = cdmath.Field("Solid temperature", cdmath.CELLS, M, 1)
  49 temperature_field_nodes = cdmath.Field("Solid temperature", cdmath.NODES, M, 1)
  50 Tin  = 330.
  51 Tout = 300.
  52
  53 for i in range(nbCells):
  54     x = M.getCell(i).x()
  55     y = M.getCell(i).y()
  56     distance = sqrt( (x - xcentre) * (x - xcentre) + (y - ycentre) * (y - ycentre) )
  57     if distance < radius:
  58         temperature_field_cells[i] = Tin
  59     else:
  60         temperature_field_cells[i] = Tout
  61
  62 temperature_field_cells.writeMED(filename, False)
  63
  64 for i in range(nbNodes):
  65     x = M.getNode(i).x()
  66     y = M.getNode(i).y()
  67     distance = sqrt( (x - xcentre) * (x - xcentre) + (y - ycentre) * (y - ycentre) )
  68     if distance < radius:
  69         temperature_field_nodes[i] = Tin
  70     else:
  71         temperature_field_nodes[i] = Tout
  72
  73 temperature_field_nodes.writeMED(filename, False)
  74
  75 # Create boundary fields
  76 ## bottom
  77 Mbottom=M.getBoundaryGroupMesh ( "Bottom" )
  78 Mbottom.writeMED(filename, False)
  79 Mbottom.writeVTK(filename)
  80 temperature_bottom_cells=cdmath.Field("Bottom temperature",cdmath.CELLS,Mbottom)
  81 for i in range(temperature_bottom_cells.getNumberOfElements()):
  82         temperature_bottom_cells[i]=Tout
  83 temperature_bottom_cells.writeMED(filename,False)
  84 temperature_bottom_cells.writeVTK(filename,True)
  85 temperature_bottom_cells.writeCSV(filename)
  86 temperature_bottom_nodes=cdmath.Field("Bottom temperature",cdmath.NODES,Mbottom)
  87 for i in range(temperature_bottom_nodes.getNumberOfElements()):
  88         temperature_bottom_nodes[i]=Tout
  89 temperature_bottom_nodes.writeMED(filename,False)
  90 temperature_bottom_nodes.writeVTK(filename,True)
  91 temperature_bottom_nodes.writeCSV(filename)
  92 ##top
  93 Mtop=M.getBoundaryGroupMesh ( "Top" )
  94 Mtop.writeMED(filename, False)
  95 Mtop.writeVTK(filename)
  96 temperature_top_cells=cdmath.Field("Top temperature",cdmath.CELLS,Mtop)
  97 for i in range(temperature_top_cells.getNumberOfElements()):
  98         temperature_top_cells[i]=Tout
  99 temperature_top_cells.writeMED(filename,False)
 100 temperature_top_cells.writeVTK(filename,True)
 101 temperature_top_cells.writeCSV(filename)
 102 temperature_top_nodes=cdmath.Field("Top temperature",cdmath.NODES,Mtop)
 103 for i in range(temperature_top_nodes.getNumberOfElements()):
 104         temperature_top_nodes[i]=Tout
 105 temperature_top_nodes.writeMED(filename,False)
 106 temperature_top_nodes.writeVTK(filename,True)
 107 temperature_top_nodes.writeCSV(filename)
 108 ##Left
 109 Mleft=M.getBoundaryGroupMesh ( "Left" )
 110 Mleft.writeMED(filename, False)
 111 Mleft.writeVTK(filename)
 112 temperature_left_cells=cdmath.Field("Left temperature",cdmath.CELLS,Mleft)
 113 for i in range(temperature_left_cells.getNumberOfElements()):
 114         temperature_left_cells[i]=Tout
 115 temperature_left_cells.writeMED(filename,False)
 116 temperature_left_cells.writeVTK(filename,True)
 117 temperature_left_cells.writeCSV(filename)
 118 temperature_left_nodes=cdmath.Field("Left temperature",cdmath.NODES,Mleft)
 119 for i in range(temperature_left_nodes.getNumberOfElements()):
 120         temperature_left_nodes[i]=Tout
 121 temperature_left_nodes.writeMED(filename,False)
 122 temperature_left_nodes.writeVTK(filename,True)
 123 temperature_left_nodes.writeCSV(filename)
 124 ##Right
 125 Mright=M.getBoundaryGroupMesh ( "Right" )
 126 Mright.writeMED(filename, False)
 127 Mright.writeVTK(filename)
 128 temperature_right_cells=cdmath.Field("Right temperature",cdmath.CELLS,Mright)
 129 for i in range(temperature_right_cells.getNumberOfElements()):
 130         temperature_right_cells[i]=Tout
 131 temperature_right_cells.writeMED(filename,False)
 132 temperature_right_cells.writeVTK(filename,True)
 133 temperature_right_cells.writeCSV(filename)
 134 temperature_right_nodes=cdmath.Field("Right temperature",cdmath.NODES,Mright)
 135 for i in range(temperature_right_nodes.getNumberOfElements()):
 136         temperature_right_nodes[i]=Tout
 137 temperature_right_nodes.writeMED(filename,False)
 138 temperature_right_nodes.writeVTK(filename,True)
 139 temperature_right_nodes.writeCSV(filename)
 140
 141 # Create fluid temperature fields
 142 temperature_field_cells = cdmath.Field("Fluid temperature", cdmath.CELLS, M, 1)
 143 temperature_field_nodes = cdmath.Field("Fluid temperature", cdmath.NODES, M, 1)
 144 Tfluid  = 300.
 145
 146 for i in range(nbCells):
 147     x = M.getCell(i).x()
 148     y = M.getCell(i).y()
 149     temperature_field_cells[i] = Tfluid
 150
 151 temperature_field_cells.writeMED(filename, False)
 152
 153 for i in range(nbNodes):
 154     x = M.getNode(i).x()
 155     y = M.getNode(i).y()
 156     temperature_field_nodes[i] = Tfluid
 157
 158 temperature_field_nodes.writeMED(filename, False)
 159
 160 # Create fluid enthalpy fields
 161 enthalpy_field_cells = cdmath.Field("Fluid enthalpy", cdmath.CELLS, M, 1)
 162 enthalpy_field_nodes = cdmath.Field("Fluid enthalpy", cdmath.NODES, M, 1)
 163 Hin  = 1.6e6
 164 Hout = 2.e6
 165
 166 for i in range(nbCells):
 167     x = M.getCell(i).x()
 168     y = M.getCell(i).y()
 169     distance = sqrt( (x - xcentre) * (x - xcentre) + (y - ycentre) * (y - ycentre) )
 170     if distance < radius:
 171         enthalpy_field_cells[i] = Hin
 172     else:
 173         enthalpy_field_cells[i] = Hout
 174
 175 enthalpy_field_cells.writeMED(filename, False)
 176
 177 for i in range(nbNodes):
 178     x = M.getNode(i).x()
 179     y = M.getNode(i).y()
 180     distance = sqrt( (x - xcentre) * (x - xcentre) + (y - ycentre) * (y - ycentre) )
 181     if distance < radius:
 182         enthalpy_field_nodes[i] = Hin
 183     else:
 184         enthalpy_field_nodes[i] = Hout
 185
 186 enthalpy_field_nodes.writeMED(filename, False)
 187
 188 # Create heat power fields
 189 heat_field_cells = cdmath.Field("Heat power", cdmath.CELLS, M, 1)
 190 heat_field_nodes = cdmath.Field("Heat power", cdmath.NODES, M, 1)
 191 phi = 1.e7
 192
 193 for i in range(nbCells):
 194     x = M.getCell(i).x()
 195     y = M.getCell(i).y()
 196     heat_field_cells[i] = phi
 197
 198 heat_field_cells.writeMED(filename, False)
 199
 200 for i in range(nbNodes):
 201     x = M.getNode(i).x()
 202     y = M.getNode(i).y()
 203     heat_field_nodes[i] = phi
 204
 205 heat_field_nodes.writeMED(filename, False)
 206
 207 # Create pressure and velocity fields for the wave equation on CELLS
 208 p0=155e5#reference pressure in a pressurised nuclear vessel
 209
 210 pressure_field = cdmath.Field("Pressure",            cdmath.CELLS, M, 1)
 211 velocity_field = cdmath.Field("Velocity",            cdmath.CELLS, M, 3)
 212
 213 for i in range(nbCells):
 214     x = M.getCell(i).x()
 215     y = M.getCell(i).y()
 216     distance = sqrt( (x - xcentre) * (x - xcentre) + (y - ycentre) * (y - ycentre) )
 217
 218     velocity_field[i,0] = 0
 219     velocity_field[i,1] = 0
 220     velocity_field[i,2] = 0
 221
 222     if distance < radius:
 223         pressure_field[i] = p0
 224         pass
 225     else:
 226         pressure_field[i] = p0/2
 227         pass
 228     pass
 229
 230 pressure_field.writeMED(filename, False)
 231 velocity_field.writeMED(filename, False)
 232
 233 #### Delete mesh and still save field in a med file
 234 m = mc.MEDCouplingCMesh()
 235 x = mc.DataArrayDouble([0.,1.,2.])
 236 y = mc.DataArrayDouble([0.,1.,2.])
 237 m.setCoords(x, y)
 238 m = m.buildUnstructured()
 239 m.setName("mesh")
 240
 241 f = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
 242 f.setMesh(m)
 243 f.setName("F")
 244
 245 da = mc.DataArrayDouble([1,2,3,4])
 246 f.setArray(da)
 247 f.setTime(0.0,0,0)
 248
 249
 250 # Maillage d'abord:
 251 fName = "/tmp/michael.med"
 252 mc.WriteUMesh(fName, m, True)
 253
 254
 255 # Maintenant juste les champs:
 256 ff = mc.MEDFileField1TS()
 257 ff.setFieldNoProfileSBT(f)
 258 ff.write(fName, 0)
 259
 260 # Tue le maillage
 261 m = 0
 262 del m
 263 import gc
 264 gc.collect()  # Make sure Python interp has called mesh destructor ...
 265
 266 # Ecrit encore du champ:
 267 da2 = da.deepCopy()
 268 f.setTime(1.0,1,0)
 269 da2 *= -1
 270 print(da2.getValues())
 271 f.setArray(da2)
 272
 273 ff = mc.MEDFileField1TS()
 274 ff.setFieldNoProfileSBT(f)  # le maillage n'existe plus, tant pis :-)
 275 ff.write(fName, 0)    # yes 0