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")
  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[0]=="Top"
  23 assert boundaryFaceGroupNames[1]=="Right"
  24 assert boundaryFaceGroupNames[2]=="Left"
  25 assert boundaryFaceGroupNames[3]=="Bottom"
  26 assert boundaryFaceGroupNames[4]=="Boundary"
  27
  28 assert boundaryNodeGroupNames[0]=="Top"
  29 assert boundaryNodeGroupNames[1]=="Right"
  30 assert boundaryNodeGroupNames[2]=="Left"
  31 assert boundaryNodeGroupNames[3]=="Bottom"
  32 assert boundaryNodeGroupNames[4]=="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 scalar fields
  48 temperature_field_cells = cdmath.Field("Temperature", cdmath.CELLS, M, 1)
  49 temperature_field_nodes = cdmath.Field("Temperature", cdmath.NODES, M, 1)
  50 Tin  = 300
  51 Tout = 400
  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 #### Delete mesh and still save field in a med file
  76 m = mc.MEDCouplingCMesh()
  77 x = mc.DataArrayDouble([0.,1.,2.])
  78 y = mc.DataArrayDouble([0.,1.,2.])
  79 m.setCoords(x, y)
  80 m = m.buildUnstructured()
  81 m.setName("mesh")
  82
  83 f = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
  84 f.setMesh(m)
  85 f.setName("F")
  86
  87 da = mc.DataArrayDouble([1,2,3,4])
  88 f.setArray(da)
  89 f.setTime(0.0,0,0)
  90
  91 # Maillage d'abord:
  92 fName = "./michael.med"
  93 mc.WriteUMesh(fName, m, True)
  94
  95 # Maintenant juste les champs:
  96 ff = mc.MEDFileField1TS()
  97 ff.setFieldNoProfileSBT(f)
  98 ff.write(fName, 0)
  99
 100 # Tue le maillage
 101 m = 0
 102 del m
 103 import gc
 104 gc.collect()  # Make sure Python interp has called mesh destructor ...
 105
 106 # Ecrit encore du champ:
 107 da2 = da.deepCopy()
 108 f.setTime(1.0,1,0)
 109 da2 += 1
 110 print(da2.getValues())
 111 f.setArray(da2)
 112
 113 ff = mc.MEDFileField1TS()
 114 ff.setFieldNoProfileSBT(f)  # le maillage n'existe plus, tant pis :-)
 115 ff.write(fName, 0)
 116 ff.write(fName, 0)  # 1 append