The most common usage is to write dedicated code (C++ or Python) linking to the library. However a graphical user interface is also available; for this please refer to the MED module documentation.
- If you don't know where to start, reading the \ref start "getting started" section and then
-taking a look at the <a class="el" href="tutorial/index.html">tutorial</a> is probably a good way to go.
+taking a look at the <a class="el" href="../tutorial/index.html">tutorial</a> is probably a good way to go.
- If you are looking for a very specific point (how can I build a mesh from scratch,
how do I write my data to a file, etc ...), taking a look at the \ref faq "FAQ" or at the numerous
\ref examples "code examples" might help.
- \ref library
- \ref functionalities
- \ref python-api
-- <a class="el" href="tutorial/index.html">Tutorial - MEDCoupling/MEDLoader in Python</a>
+- <a class="el" href="../tutorial/index.html">Tutorial - MEDCoupling/MEDLoader in Python</a>
- \ref reference
- \ref medcoupling
- \ref medloader
If you are completely new to MED, this page will help you grasp the main concepts
used overall in the \ref library "MED world", and have an idea of what you can achieve with MED.
-The <a class="el" href="tutorial/index.html">tutorial</a> is also a good way to start.
+The <a class="el" href="../tutorial/index.html">tutorial</a> is also a good way to start.
Once you are familiar with those concepts, more detailed explanations are available
in the \ref reference "reference manual".
::
- import MEDCoupling as mc
+ import medcoupling as mc
from MEDCouplingCorba import MEDCouplingUMeshServant
# Creating a mesh
import CORBA
orb = CORBA.ORB_init()
ior = orb.object_to_string(ref_m)
- print ior
+ print(ior)
# Displaying it in ParaVis
import salome
salome.salome_init()
- print "About to import module 'pvsimple' ..."
+ print("About to import module 'pvsimple' ...")
import pvsimple as pvs
- print "Module 'pvsimple' was imported!"
+ print("Module 'pvsimple' was imported!")
# From here, we use the standard ParaView API:
src1 = pvs.ParaMEDCorbaPluginSource()
nbOfCells = (N-1)*(N-1)*(N-1)
nbOfCells2D = (N-1)*(N-1)
- print "1 ********************"
+ print("1 ********************")
# Initialisation of coordinates
coordinates = []
for k in range(N):
coordinates.append(float(j))
coordinates.append(float(k))
- print "2 ********************"
+ print("2 ********************")
# Creation of meshing : need following initialisations
# => Definition of the mesh dimension
# => Definition of number of cells
mesh.allocateCells(nbOfCells+nbOfCells2D)
mesh.setName("3Dcube")
- print "3 ********************"
+ print("3 ********************")
# One remark : only one dimension cells by meshing
# Construction of volumic meshing
# => Definition of connectivity
connectivity.append(inode-N*N-N)
connectivity.append(inode-N*N-N+1)
connectivity.append(inode-N*N+1)
- print len(connectivity)
- print 8*(nbOfCells)
+ print(len(connectivity))
+ print(8*(nbOfCells))
- print "4 ********************"
+ print("4 ********************")
# Adding cells in meshing
for i in range(nbOfCells):
mesh.insertNextCell(NORM_HEXA8,8,connectivity[8*i:8*(i+1)])
pass
- print "5 ********************"
+ print("5 ********************")
# Settings of coordinates and verify if it's OK
myCoords = DataArrayDouble.New()
myCoords.setValues(coordinates,nbOfNodes,3)
mesh.setCoords(myCoords)
mesh.checkConsistencyLight()
- print "6 ********************"
+ print("6 ********************")
# Extraction of surfacic meshing
pt=[0.,0.,0.]
vec=[0.,0.,1.]
nodes = mesh.findNodesOnPlane(pt,vec,1e-12)
mesh2D = mesh.buildFacePartOfMySelfNode(nodes,True)
- #print mesh2D
+ #print(mesh2D)
mesh2D.setName("3Dcube")
mesh2D.checkConsistencyLight()
- print "7 ********************"
+ print("7 ********************")
# Creation of field : with following definition
# => Definition of the mesh support
# => Definition of field name
myCoords=DataArrayDouble.New()
sampleTab=[]
bar = mesh.computeCellCenterOfMass()
- print bar.getNbOfElems()
+ print(bar.getNbOfElems())
for i in range(nbOfCells):
x = bar.getIJ(i+1,1)
y = bar.getIJ(i+1,2)
pt=[1]
m2 = meshU.buildPartOfMySelf(pt,True);
ret,tabIdCells = meshU.areCellsIncludedIn(m2,0)
- print ret
- print tabIdCells
+ print(ret)
+ print(tabIdCells)
# Definition of the name group
tabIdCells.setName("meshGroup")
fmeshU.setName("Grid")
fmeshU.setDescription("IHopeToConvinceLastMEDMEMUsers")
myCoords = meshU.getCoords()
- print myCoords
+ print(myCoords)
fmeshU.setCoords(myCoords)
- print "**************************"
+ print("**************************")
fmeshU.setMeshAtLevel(0,meshU)
- print "**************************"
+ print("**************************")
fmeshU.setGroupsAtLevel(0,[tabIdCells],False)
- print "**************************"
+ print("**************************")
medFileName = "MEDCoupling_Gridcube3D.med"
fmeshU.write(medFileName,2)
::
- import MEDCoupling as mc
+ import medcoupling as mc
import math
# Building the coordinates of the initial hexagon, centered at 0,0
d = mc.DataArrayDouble(6,2)
d[:,0] = 3.
- d[:,1] = range(6)
+ d[:,1] = list(range(6))
d[:,1] *= math.pi/3.
d = d.fromPolarToCart()
d.setInfoOnComponents(["X [m]","Y [m]"])
- print d.getValues()
- print d
- print "Uniform array?", d.magnitude().isUniform(3.,1e-12)
+ print(d.getValues())
+ print(d)
+ print("Uniform array?", d.magnitude().isUniform(3.,1e-12))
# Translating the 7 hexagons with a translation
radius = 3.
translationToPerform = [[0.,0.],[3./2.*radius,-radius*math.sqrt(3.)/2],[3./2.*radius,radius*math.sqrt(3.)/2],[0.,radius*math.sqrt(3.)],[-3./2.*radius,radius*math.sqrt(3.)/2],[-3./2.*radius,-radius*math.sqrt(3.)/2],[0.,-radius*math.sqrt(3.)]]
oldNbOfTuples = d2.getNumberOfTuples()
c,cI = d2.findCommonTuples(1e-12)
tmp = c[cI[0]:cI[0+1]]
- print tmp
+ print(tmp)
a = cI.deltaShiftIndex()
b = a - 1
myNewNbOfTuples = oldNbOfTuples - sum(b.getValues())
o2n, newNbOfTuples = mc.DataArrayInt.ConvertIndexArrayToO2N(oldNbOfTuples,c,cI)
- print "Have I got the right number of tuples?"
- print "myNewNbOfTuples = %d, newNbOfTuples = %d" % (myNewNbOfTuples, newNbOfTuples)
+ print("Have I got the right number of tuples?")
+ print("myNewNbOfTuples = %d, newNbOfTuples = %d" % (myNewNbOfTuples, newNbOfTuples))
assert(myNewNbOfTuples == newNbOfTuples)
# Extracting the unique set of tuples
d3 = d2.renumberAndReduce(o2n, newNbOfTuples)
n2o = o2n.invertArrayO2N2N2O(newNbOfTuples)
d3_bis = d2[n2o]
- print "Are d3 and d3_bis equal ? %s" % (str(d3.isEqual(d3_bis, 1e-12)))
+ print("Are d3 and d3_bis equal ? %s" % (str(d3.isEqual(d3_bis, 1e-12))))
# Now translate everything
d3 += [3.3,4.4]
# And build an unstructured mesh representing the final pattern
m = mc.MEDCouplingUMesh("My7hexagons",2)
m.setCoords(d3)
- print "Mesh dimension is", m.getMeshDimension()
- print "Spatial dimension is", m.getCoords().getNumberOfComponents()
+ print("Mesh dimension is", m.getMeshDimension())
+ print("Spatial dimension is", m.getCoords().getNumberOfComponents())
m.allocateCells(7)
- for i in xrange(7):
+ for i in list(range(7)):
cell_connec = o2n[6*i:6*(i+1)]
m.insertNextCell(mc.NORM_POLYGON, cell_connec.getValues())
pass
::
- import MEDCoupling as mc
+ import medcoupling as mc
# Create an unstructured mesh from a Cartesian one
xarr = mc.DataArrayDouble.New(11,1)
f2.setMesh(mesh)
f2.setName("MyField2")
f2.fillFromAnalytic(1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)") # 1 means that the field should have one component
- print "Are f and f2 equal?", f.isEqualWithoutConsideringStr(f2,1e-12,1e-12)
+ print("Are f and f2 equal?", f.isEqualWithoutConsideringStr(f2,1e-12,1e-12))
#
da1 = f.getArray() # a DataArrayDouble, which is a direct reference (not a copy) of the field's values
ids1 = da1.findIdsInRange(0.,5.)
# Check that fPart1Cpy and fPart1 are the same
fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
fPart1Cpy.getArray().abs()
- print "Are the fields equal?", (fPart1Cpy.getArray().accumulate()[0]<1e-12)
+ print("Are the fields equal?", (fPart1Cpy.getArray().accumulate()[0]<1e-12))
# Aggregate fields
fPart12 = mc.MEDCouplingFieldDouble.MergeFields([fPart1,fPart2])
fPart12.writeVTK("ExoField_fPart12.vtu")
arr2 = f.getValueOnMulti(bary)
delta = arr1-arr2
delta.abs()
- print "Is field evaluation matching?", (delta.accumulate()[0]<1e-12)
+ print("Is field evaluation matching?", (delta.accumulate()[0]<1e-12))
# ExtensiveMaximum computations
integ1 = fPart12.integral(0,True)
integ1_bis = fPart12.getArray().accumulate()[0]
- print "First integral matching ?", ( abs(integ1 - integ1_bis) < 1e-8 )
+ print("First integral matching ?", ( abs(integ1 - integ1_bis) < 1e-8 ))
fPart12.getMesh().scale([0.,0.,0.], 1.2)
integ2 = fPart12.integral(0,True)
- print "Second integral matching ?", ( abs(integ2-integ1_bis*1.2*1.2*1.2) < 1e-8 )
+ print("Second integral matching ?", ( abs(integ2-integ1_bis*1.2*1.2*1.2) < 1e-8 ))
# Explosion of field
fVec = mesh.fillFromAnalytic(mc.ON_CELLS,3,"(x-5.)*IVec+(y-5.)*JVec+(z-5.)*KVec")
fVecPart1 = fVec.buildSubPart(ids1)
# Get available time steps
data = ml.MEDFileData("agitateur.med")
ts = data.getFields()[0].getTimeSteps()
- print ts
+ print(ts)
# Get position of the swirler
fMts = data.getFields()["DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM"]
f1ts = fMts[(2,-1)]
torquePerCellOnSkin = ml.DataArrayDouble.CrossProduct(posSkin,forceVectSkin)
zeTorque = torquePerCellOnSkin.accumulate()
- print "couple = %r N.m" % zeTorque[2]
+ print("couple = %r N.m" % zeTorque[2])
# Power computation
speedMts = data.getFields()["VITESSE_ELEM_DOM"]
speed1ts = speedMts[(2,-1)]
speedOnSkin = speedMc.getArray()[tupleIdsInField]
powerSkin = ml.DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
power = powerSkin.accumulate()[0]
- print "power = %r W"%(power)
+ print("power = %r W"%(power))
# Eigen vector computation
x2 = posSkin[:,0]*posSkin[:,0]
x2 = x2.accumulate()[0]
inertiaSkinValues, inertiaSkinVects = np.linalg.eig(inertiaSkin)
pos = max(enumerate(inertiaSkinValues), key=lambda x: x[1])[0]
vect0 = inertiaSkinVects[pos].tolist()[0]
- print vect0
+ print(vect0)
def computeAngle(locAgitateur1ts):
fMc = locAgitateur1ts.getFieldAtLevel(ml.ON_CELLS,0)
from math import acos, sqrt
angle2 = len(ts)*[0.]
- for pos in xrange(2,len(vects)):
+ for pos in list(range(2,len(vects))):
norm1 = sqrt(vects[pos-1][0]*vects[pos-1][0]+vects[pos-1][1]*vects[pos-1][1])
norm2 = sqrt(vects[pos][0]*vects[pos][0]+vects[pos][1]*vects[pos][1])
crs = vects[pos-1][0]*vects[pos][0]+vects[pos-1][1]*vects[pos][1]
pass
omega=sum(angle2)/(ts[-1][2]-ts[0][2])
- print sum(angle2)
+ print(sum(angle2))
- print "At timestep (%d,%d) (physical time=%r s) the torque is: %r N.m, power/omega=%r N.m " % (ts[2][0],ts[2][1],ts[2][2],zeTorque[2],power/omega)
+ print("At timestep (%d,%d) (physical time=%r s) the torque is: %r N.m, power/omega=%r N.m " % (ts[2][0],ts[2][1],ts[2][2],zeTorque[2],power/omega))
# Read and clean Fixe.med
fixe = ml.MEDFileMesh.New("Fixe.med")
fixm = fixe.getMeshAtLevel(0)
- print "Nb of nodes in the file : %i " % (fixm.getNumberOfNodes())
+ print("Nb of nodes in the file : %i " % (fixm.getNumberOfNodes()))
fixm.mergeNodes(1e-10)
- print "Nb of non duplicated nodes : %i" % (fixm.getNumberOfNodes())
+ print("Nb of non duplicated nodes : %i" % (fixm.getNumberOfNodes()))
# Read and clean Mobile.med
mobile = ml.MEDFileMesh.New("Mobile.med")
mobm = mobile.getMeshAtLevel(0)
mobm2.writeVTK("mobm2.vtu")
# mobm2 is in several pieces, take the first one
zonesInMobm = mobm.partitionBySpreadZone()
- print "Nb of zones in mobm : %i" % (len(zonesInMobm))
+ print("Nb of zones in mobm : %i" % (len(zonesInMobm)))
zone1Mobm = mobm[zonesInMobm[0]]
zone1Mobm.zipCoords()
displayVTK(zone1Mobm, "zone1Mobm.vtu")
displayVTK(partFixmWithoutZone1Mobm,"partFixmWithoutZone1Mobm.vtu")
# Check that intersection worked properly
# Check #0
- areaPartFixm = partFixm.getMeasureField(ml.ON_CELLS).getArray()
+ areaPartFixm = partFixm.getMeasureField(False).getArray() # if set to True returns the absolut field value
areaPartFixm.abs()
- areaPartFixMob = partFixMob.getMeasureField(ml.ON_CELLS).getArray()
+ areaPartFixMob = partFixMob.getMeasureField(False).getArray()
areaPartFixMob.abs()
val1=areaPartFixm.accumulate()[0]
val2=areaPartFixMob.accumulate()[0]
- print "Check #0 %lf == %lf with precision 1e-8? %s" % (val1,val2,str(abs(val1-val2)<1e-8))
+ print("Check #0 %lf == %lf with precision 1e-8? %s" % (val1,val2,str(abs(val1-val2)<1e-8)))
# Check #1
- areaZone1Mobm = zone1Mobm.getMeasureField(ml.ON_CELLS).getArray()
+ areaZone1Mobm = zone1Mobm.getMeasureField(False).getArray()
areaZone1Mobm.abs()
val3 = areaZone1Mobm.accumulate()[0]
ids4 = iMob.findIdsNotEqual(-1)
areaPartFixMob2 = areaPartFixMob[ids4]
val4 = areaPartFixMob2.accumulate()[0]
- print "Check #1 %lf == %lf with precision 1e-8 ? %s" % (val3,val4,str(abs(val3-val4)<1e-8))
+ print("Check #1 %lf == %lf with precision 1e-8 ? %s" % (val3,val4,str(abs(val3-val4)<1e-8)))
# Check #2
isCheck2OK = True
- for icell in xrange(partFixm.getNumberOfCells()):
+ for icell in list(range(partFixm.getNumberOfCells())):
ids5 = iPart.findIdsEqual(icell)
areaOfCells = areaPartFixMob[ids5]
areaOfCells.abs()
isCheck2OK = False
pass
pass
- print "Check #2? %s" % (str(isCheck2OK))
+ print("Check #2? %s" % (str(isCheck2OK)))
# Indicator field creation
f = ml.MEDCouplingFieldDouble(ml.ON_CELLS,ml.ONE_TIME)
m = partFixMob.deepCopy()
::
- import MEDCoupling as mc
+ import medcoupling as mc
#
# NumPy
arr = mc.DataArrayDouble(12)
arr[:] = 4.
nparr = arr.toNumPyArray()
- print nparr.__repr__()
- print nparr.tolist()
+ print(nparr.__repr__())
+ print(nparr.tolist())
nparr[::2] = 7.
- print nparr.__repr__()
- print arr.__repr__()
+ print(nparr.__repr__())
+ print(arr.__repr__())
del arr
import gc; gc.collect() # Make sure the object has been deleted
- print nparr.__repr__()
+ print(nparr.__repr__())
arr2 = mc.DataArrayDouble(nparr)
- print arr2.__repr__()
+ print(arr2.__repr__())
nparr[:] = 5.
- print nparr.__repr__()
- print arr2.__repr__()
+ print(nparr.__repr__())
+ print(arr2.__repr__())
# Writing to file
f = open("toto.data","w+b")
a = np.memmap(f,dtype='float64',mode='w+',offset=0,shape=nparr.shape)
a[:] = 3.14
f.flush()
b = np.memmap(f2,dtype='float64',mode='r',offset=0,shape=(12,))
- print b.__repr__()
+ print(b.__repr__())
#
# SciPy
#
numberOfNodes = 25
numberOfCells = 12
- print "1 ********************"
+ print("1 ********************")
spaceDimension = 2
# Coordinates of central polygon
coordinates.append(xtmp+X[(i+1)%6])
coordinates.append(ytmp+Y[(i+1)%6])
- print "2 ********************"
+ print("2 ********************")
# Creation of mesh
mesh=MEDCouplingUMesh.New()
mesh.setMeshDimension(2)
myCoords.setValues(coordinates,numberOfNodes,2)
mesh.setCoords(myCoords)
- print "3 ********************"
+ print("3 ********************")
# Connectivity of triangular meshing
connectivity = []
for i in range(6):
mesh.insertNextCell(NORM_TRI3,3,connectivity[3*i:3*(i+1)])
pass
- print "4 ********************"
+ print("4 ********************")
# Connectivity of hexagons
connectivity = []
for i in range(6):
mesh.insertNextCell(NORM_POLYGON,6,connectivity[6*i:6*(i+1)])
pass
- print "5 ********************"
+ print("5 ********************")
mesh.checkConsistencyLight()
medFileName = "MEDCoupling_Fleur.med"
sinus = sin(d)
if abs(res[0]-sinus)<1.e-5:
- print "OK"
+ print("OK")
else:
- print "KO"
+ print("KO")
::
- import MEDCoupling as mc
+ import medcoupling as mc
from MEDCouplingRemapper import MEDCouplingRemapper
# Target mesh
arr = mc.DataArrayDouble(11)
remap.prepare(srcMesh,trgMesh,"P0P0")
# Check matrix
myMatrix = remap.getCrudeMatrix()
- print myMatrix
+ print(myMatrix)
sumByRows = mc.DataArrayDouble(len(myMatrix))
for i,wIt in enumerate(sumByRows):
su = 0.
for it in myMatrix[i]:
su += myMatrix[i][it]
wIt[0] = su
- print "Is interpolation well prepared?", sumByRows.isUniform(1.,1e-12)
+ print("Is interpolation well prepared?", sumByRows.isUniform(1.,1e-12))
# Source field construction
srcField = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
srcField.setMesh(srcMesh)
# IntensiveMaximum
integSource = srcField.integral(True)[0]
integTarget = trgFieldCV.integral(True)[0]
- print "IntensiveMaximum -- integrals: %lf == %lf" % (integSource, integTarget)
+ print("IntensiveMaximum -- integrals: %lf == %lf" % (integSource, integTarget))
accSource = srcField.getArray().accumulate()[0]
accTarget = trgFieldCV.getArray().accumulate()[0]
- print "IntensiveMaximum -- sums: %lf != %lf" % (accSource, accTarget)
+ print("IntensiveMaximum -- sums: %lf != %lf" % (accSource, accTarget))
# ExtensiveConservation
srcField.setNature(mc.ExtensiveConservation)
trgFieldI = remap.transferField(srcField,1e300)
#
integSource = srcField.integral(True)[0]
integTarget = trgFieldI.integral(True)[0]
- print "ExtensiveConservation -- integrals: %lf != %lf" % (integSource, integTarget)
-
+ print("ExtensiveConservation -- integrals: %lf != %lf" % (integSource, integTarget))
accSource = srcField.getArray().accumulate()[0]
accTarget = trgFieldI.getArray().accumulate()[0]
- print "ExtensiveConservation -- sums: %lf == %lf" % (accSource, accTarget)
+ print("ExtensiveConservation -- sums: %lf == %lf" % (accSource, accTarget))
::
- import MEDCoupling as mc
+ import medcoupling as mc
# Build a 3D mesh from scratch mixing HEXA8 and POLYHED
coords=[0.,0.,0., 1.,1.,0., 1.,1.25,0., 1.,0.,0., 1.,1.5,0., 2.,0.,0., 2.,1.,0., 1.,2.,0., 0.,2.,0., 3.,1.,0.,
mesh3D.orientCorrectlyPolyhedrons()
mesh3D.sortCellsInMEDFileFrmt()
mesh3D.checkConsistencyLight()
- renum = mc.DataArrayInt(60) ; renum[:15]=range(15,30) ; renum[15:30]=range(15) ; renum[30:45]=range(45,60) ; renum[45:]=range(30,45)
+ renum = mc.DataArrayInt(60) ; renum[:15]=list(range(15,30)) ; renum[15:30]=list(range(15)) ; renum[30:45]=list(range(45,60)) ; renum[45:]=list(range(30,45))
mesh3D.renumberNodes(renum,60)
# Scale coordinates from meters to centimeters
mesh3D.getCoords()[:] *= 100.
n_cells = mesh2D.getNumberOfCells()
cellIdsSol3 = extMesh.getMesh3DIds()[n_cells:2*n_cells]
# Compare the 3 methods
- print cellIdsSol1.getValues()
- print cellIdsSol2.getValues()
- print cellIdsSol3.getValues()
+ print(cellIdsSol1.getValues())
+ print(cellIdsSol2.getValues())
+ print(cellIdsSol3.getValues())
# Extract part of the mesh
mesh3DPart = mesh3D[cellIdsSol2] # equivalent to mesh3DPart = mesh3D.buildPartOfMySelf(cellIdsSol2,True)
mesh3DPart.zipCoords()
# Check geometric type ordering
- #print mesh3DPart.advancedRepr()
- print mesh3DPart.checkConsecutiveCellTypesAndOrder([mc.NORM_HEXA8,mc.NORM_POLYHED])
- print mesh3DPart.checkConsecutiveCellTypes()
- #print mesh3DPart.advancedRepr()
+ #print(mesh3DPart.advancedRepr())
+ print(mesh3DPart.checkConsecutiveCellTypesAndOrder([mc.NORM_HEXA8,mc.NORM_POLYHED]))
+ print(mesh3DPart.checkConsecutiveCellTypes())
+ #print(mesh3DPart.advancedRepr())
# Extract cells along a line - Solution 1
baryXY = bary[:,[0,1]]
baryXY -= [250.,150.]
magn = bary2.magnitude()
ids = magn.findIdsInRange(0.,1e-12)
idStart = int(ids) # ids is assumed to contain only one value, if not an exception is thrown
- ze_range = range(idStart,mesh3D.getNumberOfCells(),mesh2D.getNumberOfCells())
+ ze_range = list(range(idStart,mesh3D.getNumberOfCells(),mesh2D.getNumberOfCells()))
cellIds2Sol2 = extMesh.getMesh3DIds()[ze_range]
# Construct the final sub-part
mesh3DSlice2 = mesh3D[cellIds2Sol1]
meshMEDFileRead = ml.MEDFileMesh.New("TargetMesh2.med") # a new is needed because it returns a MEDFileUMesh (MEDFileMesh is abstract)
meshRead0 = meshMEDFileRead.getMeshAtLevel(0)
meshRead1 = meshMEDFileRead.getMeshAtLevel(-1)
- print "Is level 0 in the file equal to 'targetMesh'?", meshRead0.isEqual(targetMesh,1e-12)
- print "Is level 0 in the file equal to 'targetMesh1'?", meshRead1.isEqual(targetMesh1,1e-12)
+ print("Is level 0 in the file equal to 'targetMesh'?", meshRead0.isEqual(targetMesh,1e-12))
+ print("Is level 0 in the file equal to 'targetMesh1'?", meshRead1.isEqual(targetMesh1,1e-12))
# Read groups
- print meshMEDFileRead.getGrpNonEmptyLevels("grp0_Lev0")
+ print(meshMEDFileRead.getGrpNonEmptyLevels("grp0_Lev0"))
grp0_0_read = meshMEDFileRead.getGroupArr(0,"grp0_Lev0")
- print "Is group 'grp0_Lev0' equal to what is read in the file?" , grp0_0_read.isEqual(grp0_0)
+ print("Is group 'grp0_Lev0' equal to what is read in the file?" , grp0_0_read.isEqual(grp0_0))
#
# Fields
#
fMEDFileRead = ml.MEDFileField1TS("TargetMesh2.med",f.getName(),7,8)
fRead1 = fMEDFileRead.getFieldOnMeshAtLevel(ml.ON_CELLS,0,meshMEDFileRead) # Quickest way, not re-reading mesh in the file.
fRead2 = fMEDFileRead.getFieldAtLevel(ml.ON_CELLS,0) # Like above, but this time the mesh is read!
- print "Does the field remain OK with the quick method?", fRead1.isEqual(f,1e-12,1e-12)
- print "Does the field remain OK with the slow method?", fRead2.isEqual(f,1e-12,1e-12)
+ print("Does the field remain OK with the quick method?", fRead1.isEqual(f,1e-12,1e-12))
+ print("Does the field remain OK with the slow method?", fRead2.isEqual(f,1e-12,1e-12))
#
# Writing and Reading fields on profile using MEDLoader advanced API
#
#
fMEDFileRead2 = ml.MEDFileField1TS("TargetMesh2.med",fPart.getName(),7,8)
fPartRead, pflRead = fMEDFileRead2.getFieldWithProfile(ml.ON_CELLS,0,meshMEDFileRead)
- print "Is the partial field correctly read?", fPartRead.isEqualWithoutConsideringStr(fPart.getArray(),1e-12)
- print "Is the list of cell identifiers matching?", pflRead.isEqualWithoutConsideringStr(pfl)
+ print("Is the partial field correctly read?", fPartRead.isEqualWithoutConsideringStr(fPart.getArray(),1e-12))
+ print("Is the list of cell identifiers matching?", pflRead.isEqualWithoutConsideringStr(pfl))
ml.WriteUMesh("TargetMesh.med",targetMesh,True) # True means 'from scratch'
# Re-read it and test equality
meshRead = ml.ReadUMeshFromFile("TargetMesh.med",targetMesh.getName(),0)
- print "Is the read mesh equal to 'targetMesh' ?", meshRead.isEqual(targetMesh,1e-12)
+ print("Is the read mesh equal to 'targetMesh' ?", meshRead.isEqual(targetMesh,1e-12))
# Writing a field and its support mesh in one go
f = ml.MEDCouplingFieldDouble.New(ml.ON_CELLS, ml.ONE_TIME)
f.setTime(5.6,7,8) # Declare the timestep associated to the field
ml.WriteField("MyFirstField.med",f,True)
# Re-read it and test equality
f2 = ml.ReadFieldCell("MyFirstField.med", f.getMesh().getName(), 0, f.getName(), 7, 8)
- print "Is the read field identical to 'f' ?", f2.isEqual(f,1e-12,1e-12)
+ print("Is the read field identical to 'f' ?", f2.isEqual(f,1e-12,1e-12))
# Writing in several steps
ml.WriteUMesh("MySecondField.med",f.getMesh(),True)
ml.WriteFieldUsingAlreadyWrittenMesh("MySecondField.med",f)
ml.WriteFieldUsingAlreadyWrittenMesh("MySecondField.med",f2)
# Re-read and test this two-timestep field
f3 = ml.ReadFieldCell("MySecondField.med",f.getMesh().getName(),0,f.getName(),7,8)
- print "Is the field read in file equals to 'f' ?", f.isEqual(f3,1e-12,1e-12)
+ print("Is the field read in file equals to 'f' ?", f.isEqual(f3,1e-12,1e-12))
f4 = ml.ReadFieldCell("MySecondField.med",f.getMesh().getName(),0,f.getName(),9,10)
- print "Is the field read in file equals to 'f2' ?", f2.isEqual(f4,1e-12,1e-12)
+ print("Is the field read in file equals to 'f2' ?", f2.isEqual(f4,1e-12,1e-12))
cellFieldCpy = cellField.deepCopy()
cellFieldCpy.substractInPlaceDM(cellField_read,10,1e-12)
cellFieldCpy.getArray().abs()
- print cellFieldCpy.getArray().isUniform(0.,1e-12)
+ print(cellFieldCpy.getArray().isUniform(0.,1e-12))
#
nodeField0_read = ml.ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
nodeField1_read = ml.ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
nodeFieldCpy = nodeField.deepCopy()
nodeFieldCpy.mergeNodes(1e-10)
nodeFieldCpy.substractInPlaceDM(nodeField_read,10,1e-12)
- print nodeFieldCpy.getArray().isUniform(0.,1e-12)
+ print(nodeFieldCpy.getArray().isUniform(0.,1e-12))
#
# Merging - Optimal method
#
myCoords=DataArrayDouble.New()
sampleTab=[]
bar = mesh.computeCellCenterOfMass()
- print bar.getNbOfElems()
+ print(bar.getNbOfElems())
for i in range(nbOfCells):
x = bar.getIJ(...)
y = bar.getIJ(...)
vec=[0.,0.,1.]
nodes = mesh.findNodesOnPlane(pt,vec,1e-12)
mesh2D = mesh.buildFacePartOfMySelfNode(nodes,True)
- #print mesh2D
+ #print(mesh2D)
mesh2D.setName("3Dcube")
mesh2D.checkConsistencyLight()
sinus = sin(d)
if abs(res[0]-sinus)<1.e-5:
- print "OK"
+ print("OK")
else:
- print "KO"
+ print("KO")
Solution
~~~~~~~~
Starting the implementation
~~~~~~~~~~~~~~~~~~~~~~~~~~~
-To start with, import the whole python module MEDCoupling and the standard math module. ::
+To start with, import the whole python module medcoupling and the standard math module. ::
- from MEDCoupling import *
+ from medcoupling import *
import math
This makes the following available:
Initialize the 2nd component of each tuple i with the index value i. ::
- d[:,1]=range(6)
+ d[:,1]=list(range(6))
Multiply the 2nd component by pi/3. ::
Display the values only in form of a Python list. ::
- print d.getValues()
+ print(d.getValues())
Display d. ::
- print d
+ print(d)
Verify that for each tuple of d its norm (returned by the magnitude() method) is indeed 3 (with an error less than 1.e-12):
::
- print d.magnitude().isUniform(3.,1e-12)
+ print(d.magnitude().isUniform(3.,1e-12))
Duplication and aggregation of DataArrayDouble
An alternative (and more compact) way to do it : ::
- ds=[d.deepCopy() for i in xrange(len(translationToPerform))]
+ ds=[d.deepCopy() for i in list(range(len(translationToPerform)))]
for (elt,t) in zip(ds,translationToPerform) : elt+=t
Aggregating DataArrayDouble
Get the list of tuple ids forming the first group of common nodes and store the result in tmp. ::
tmp=c[cI[0]:cI[0+1]]
- print tmp
+ print(tmp)
Check the result: all the tuples stored in tmp point to identical coordinates in d2.::
- print d2[tmp]
+ print(d2[tmp])
.. note:: we see the tuple (3.,0.) repeated 3 times (with an error margin below 1e-12).
We get for free the number of elements in Y, i.e. the variable newNbOfTuples. ::
o2n,newNbOfTuples=DataArrayInt.ConvertIndexArrayToO2N(oldNbOfTuples,c,cI)
- print "Have I got the right result? %s"%(str(myNewNbOfTuples==newNbOfTuples))
+ print("Have I got the right result? %s"%(str(myNewNbOfTuples==newNbOfTuples)))
Using o2n and newNbOfTuples invoke DataArrayDouble.renumberAndReduce() on d2. ::
Using n2o we can deduce d3_bis from d2. ::
d3_bis=d2[n2o]
- print "Have I got the right result (2)? %s"%(str(d3.isEqual(d3_bis,1e-12)))
+ print("Have I got the right result (2)? %s"%(str(d3.isEqual(d3_bis,1e-12))))
Translate all tuples at once
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
m.allocateCells(7)
-Finally thanks to o2n we know the connectivity of all 7 hexagons using the coordinates stored in d3.
+Finally thanks to o2n we know the connectivity of all 7 hexagons using the coordinates stored in d3. ::
- for i in xrange(7):
- m.insertNextCell(NORM_POLYGON,o2n[6*i:6*(i+1)].getValues())
- pass
+ for i in list(range(7)):
+ m.insertNextCell(NORM_POLYGON,o2n[6*i:6*(i+1)].getValues())
+ pass
Check that m is coherent. ::
Début de l'implémentation
~~~~~~~~~~~~~~~~~~~~~~~~~
-Pour commencer l'exercice importer le module Python ``MEDCoupling`` et l'aliaser avec ``mc`` (ca nous évitera des noms trop longs). Importer aussi le module ``math``. ::
+Pour commencer l'exercice importer le module Python ``medcoupling`` et l'aliaser avec ``mc`` (ca nous évitera des noms trop longs). Importer aussi le module ``math``. ::
- import MEDCoupling as mc
+ import medcoupling as mc
import math
On rappelle que toutes les méthodes statiques du module commencent par une majuscule.
n'a qu'une seule composante. ::
d_example = mc.DataArrayDouble([0.0,1.0,2.5])
- print d_example
+ print(d_example )
.. note:: Le tableau ``d`` contient maintenant 12 valeurs groupées en 6 tuples contenant chacun 2 composantes.
Les valeurs dans ``d`` ne sont pas encore assignées.
Initialiser la 2ème composante de chaque tuple i avec la valeur i. ::
- d[:,1] = range(6)
+ d[:,1] = list(range(6))
Multiplier la seconde composante de chacun des tuples par pi/3. ::
Afficher ``d`` tel quel. ::
- print d
+ print(d)
Afficher juste les valeurs sous forme d'une liste python. ::
- print d.getValues()
+ print(d.getValues())
Vérifier que pour chaque tuple désormais dans ``d``, sa norme (méthode ``magnitude()``) est bien égale à 3.0, à 1.e-12 près (méthode ``isUniform()``) ::
- print "Uniform array?", d.magnitude().isUniform(3.,1e-12)
+ print("Uniform array?", d.magnitude().isUniform(3.,1e-12))
Duplication et agrégation
Une autre façon de faire un peu plus compacte (pour les amoureux des *one-liner*) : ::
- ds = [d + translationToPerform[i] for i in xrange(len(translationToPerform))]
+ ds = [d + translationToPerform[i] for i in list(range(len(translationToPerform)))]
Agrégation de tableaux
~~~~~~~~~~~~~~~~~~~~~~
Afficher ``tmp``. ::
tmp = c[cI[0]:cI[0+1]]
- print tmp
+ print(tmp)
Vérifier, en l'affichant, que pour tous les identifiants de tuples dans ``tmp``, leurs tuples sont bien égaux dans ``d2``. ::
- print d2[tmp]
+ print(d2[tmp])
.. note:: On voit que le tuple (3.,0.) à 1e-12 près est répété 3 fois et ``tmp`` donne les positions respectives de
ces 3 répétitions.
On récupère au passage card(Y) c'est-à-dire le ``newNbOfTuples``. ::
o2n, newNbOfTuples = mc.DataArrayInt.ConvertIndexArrayToO2N(oldNbOfTuples,c,cI)
- print "Have I got the right number of tuples?"
- print "myNewNbOfTuples = %d, newNbOfTuples = %d" % (myNewNbOfTuples, newNbOfTuples)
+ print("Have I got the right number of tuples?")
+ print("myNewNbOfTuples = %d, newNbOfTuples = %d" % (myNewNbOfTuples, newNbOfTuples))
assert(myNewNbOfTuples == newNbOfTuples)
Nous pouvons maintenant constuire le tableau de points uniques ``d3``. A l'aide de ``o2n``
A l'aide de ``n2o`` on peut construire un ``d3_bis`` à partir de ``d2``, et qui contient la même chose que le ``d3`` précédent. ::
d3_bis = d2[n2o]
- print "Are d3 and d3_bis equal ? %s" % (str(d3.isEqual(d3_bis, 1e-12)))
+ print("Are d3 and d3_bis equal ? %s" % (str(d3.isEqual(d3_bis, 1e-12))))
Translater tous les tuples
~~~~~~~~~~~~~~~~~~~~~~~~~~
m = mc.MEDCouplingUMesh("My7hexagons",2)
m.setCoords(d3)
- print "Mesh dimension is", m.getMeshDimension()
- print "Spatial dimension is", m.getCoords().getNumberOfComponents()
+ print("Mesh dimension is", m.getMeshDimension())
+ print("Spatial dimension is", m.getCoords().getNumberOfComponents())
Maintenant, allouer le nombre de cellules avec (un majorant du) nombre attendu de cellules. ::
Enfin grâce à ``o2n`` on a la *connectivité* (i.e. la liste des points formant un hexagone)
des 7 hexagones utilisant les coordonnées ``d3``. ::
- for i in xrange(7):
+ for i in range(7):
cell_connec = o2n[6*i:6*(i+1)]
m.insertNextCell(mc.NORM_POLYGON, cell_connec.getValues())
pass
Implementation start
~~~~~~~~~~~~~~~~~~~~
-Import the MEDCoupling Python module. ::
+Import the medcoupling Python module. ::
- from MEDCoupling import *
+ from medcoupling import *
We are going to create a MEDCouplingUMesh from a 3D cartesian mesh. Each direction will contain 10 cells and 11 nodes. The generated MEDCouplingUMesh
will contain 1000 cells. ::
Compare the two fields:
Compare f and f2 with a precision of 1e-12 on coordinates and 1e-12 on values. ::
- print "f and f2 are equal: %s"%(f.isEqualWithoutConsideringStr(f2,1e-12,1e-12))
+ print("f and f2 are equal: %s"%(f.isEqualWithoutConsideringStr(f2,1e-12,1e-12)))
Builing of a subpart of a field
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
fPart1Cpy.getArray().abs()
- print "Fields are the same? %s"%(fPart1Cpy.getArray().accumulate()[0]<1e-12)
+ print("Fields are the same? %s"%(fPart1Cpy.getArray().accumulate()[0]<1e-12))
.. note:: This is in fact a very special case of interpolation. Except that here
we assume that the supports of "fPart1" and "fPart1Cpy" are equal, discarding any
arr2=f.getValueOnMulti(bary)
delta=arr1-arr2
delta.abs()
- print "Check OK: %s"%(delta.accumulate()[0]<1e-12)
+ print("Check OK: %s"%(delta.accumulate()[0]<1e-12))
.. note:: In this context and for example for a field on cells, "evaluate" at a point means returning the value of the cell containing the point.
.. note:: This technique can be used to quickly assess the quality of an interpolation.
Début de l'implementation
~~~~~~~~~~~~~~~~~~~~~~~~~
-Importer le module Python ``MEDCoupling``. ::
+Importer le module Python ``medcoupling``. ::
- import MEDCoupling as mc
+ import medcoupling as mc
Créer un ``MEDCouplingUMesh`` à partir d'un maillage 3D cartésien. Chaque direction contiendra 10 cells
et 11 nodes. Le ``MEDCouplingUMesh`` résultant contiendra ainsi 1000 cells. ::
Comparer les deux champs : comparer ``f`` et ``f2`` avec une précision de 1e-12 sur les coordonnées et
de 1e-13 sur les valeurs. ::
- print "Are f and f2 equal?", f.isEqualWithoutConsideringStr(f2,1e-12,1e-13)
+ print("Are f and f2 equal?", f.isEqualWithoutConsideringStr(f2,1e-12,1e-13))
.. note:: Le ``WithoutConsideringStr`` dans le nom de la méthode précédente indique que les noms des champs ne seront
fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
fPart1Cpy.getArray().abs()
- print "Equal field ? %s" % (fPart1Cpy.getArray().accumulate()[0]<1e-12)
+ print("Equal field ? %s" % (fPart1Cpy.getArray().accumulate()[0]<1e-12))
.. note:: La renumérotation effectuée ici représente en fait d'un cas très particulier
d'interpolation. Effectivement l'hypothèse est faite que les supports
arr2 = f.getValueOnMulti(bary)
delta = arr1-arr2
delta.abs()
- print "Is field evaluation matching?", (delta.accumulate()[0]<1e-12)
+ print("Is field evaluation matching?", (delta.accumulate()[0]<1e-12))
.. note:: Dans ce contexte, et pour un champ aux cellules (P0) par exemple, "évaluer" en un point signifie retourner la valeur
de la cellule contenant le point donné.
integ1 = fPart12.integral(0,True)
integ1_bis = fPart12.getArray().accumulate()[0]
- print "First integral matching ?", ( abs(integ1 - integ1_bis) < 1e-8 )
+ print("First integral matching ?", ( abs(integ1 - integ1_bis) < 1e-8 ))
Ensuite appliquer une homotétie de facteur 1.2 centrée en [0.,0.,0.] sur le support de ``fPart12`` (c'est-à-dire son maillage).
Quelle est alors la nouvelle valeur de l'intégrale ? ::
fPart12.getMesh().scale([0.,0.,0.], 1.2)
integ2 = fPart12.integral(0,True)
- print "Second integral matching ?", ( abs(integ2-integ1_bis*1.2*1.2*1.2) < 1e-8 )
+ print("Second integral matching ?", ( abs(integ2-integ1_bis*1.2*1.2*1.2) < 1e-8 ))
Exploser un champ - Vecteurs de déplacement
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Implementation start
~~~~~~~~~~~~~~~~~~~~
-Import the MEDCoupling Python module. ::
+Import the medcoupling Python module. ::
- from MEDCoupling import *
+ from medcoupling import *
We now build a mesh containing artificially two types of cell (NORM_HEXA8 and NORM_POLYHED) to highlight the possibility to work with non-homogeneous cell types.
mesh3D is an extruded mesh containing 18 cells composed into 3 levels along Z of 6 cells.
mesh3D.orientCorrectlyPolyhedrons()
mesh3D.sortCellsInMEDFileFrmt()
mesh3D.checkConsistencyLight()
- renum=DataArrayInt.New(60) ; renum[:15]=range(15,30) ; renum[15:30]=range(15) ; renum[30:45]=range(45,60) ; renum[45:]=range(30,45)
+ renum=DataArrayInt.New(60)
+ renum[:15] = list(range(15,30))
+ renum[15:30] = list(range(15))
+ renum[30:45] = list(range(45,60))
+ renum[45:] = list(range(30,45))
mesh3D.renumberNodes(renum,60)
Convert coordinate unit from meters to centimeters
It is now possible to check that the 3 solutions are the same : ::
- for i in xrange(3):
- exec("print cellIdsSol%s.getValues()"%(i+1))
+ for i in list(range(3)):
+ exec("print(cellIdsSol%s.getValues())"%(i+1))
Extract a sub-part of mesh3D
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
At this point mesh3DPart only contains 30 nodes and 6 cells. To prepare to MED file I/O we have to check if mesh3DPart is ready to be written safely into a MED file (i.e. if the cells are indeed ordered by type). ::
- print mesh3DPart.checkConsecutiveCellTypesAndOrder([NORM_HEXA8,NORM_POLYHED])
+ print(mesh3DPart.checkConsecutiveCellTypesAndOrder([NORM_HEXA8,NORM_POLYHED]))
Or: ::
- print mesh3DPart.checkConsecutiveCellTypes()
+ print(mesh3DPart.checkConsecutiveCellTypes())
You can also print the content of the mesh "mesh3Dpart": ::
- print mesh3DPart.advancedRepr()
+ print(mesh3DPart.advancedRepr())
We see that mesh3DPart contains 6 cells, 4 HEXA8 then 2 POLYHED. Everything's OK: the cells are grouped by geometrical type.
magn=bary2.magnitude()
ids=magn.findIdsInRange(0.,1e-12)
idStart=int(ids) # ids is assumed to contain only one value, if not an exception is thrown
- cellIds2Sol2=extMesh.getMesh3DIds()[range(idStart,mesh3D.getNumberOfCells(),mesh2D.getNumberOfCells())]
+ cellIds2Sol2=extMesh.getMesh3DIds()[list(range(idStart,mesh3D.getNumberOfCells(),mesh2D.getNumberOfCells()))]
Now, build the sub-part of mesh3D using cell IDs in cellIds2Sol1. ::
Début de l'implémentation
~~~~~~~~~~~~~~~~~~~~~~~~~
-Importer le module Python ``MEDCoupling``. ::
+Importer le module Python ``medcoupling``. ::
- import MEDCoupling as mc
+ import medcoupling as mc
Construire un maillage. Ce maillage ``mesh3D`` contient artificiellement 2 types de cellules (``mc.NORM_HEXA8`` et ``mc.NORM_POLYHED``)
pour appréhender le mélange de types geometriques.
mesh3D.orientCorrectlyPolyhedrons()
mesh3D.sortCellsInMEDFileFrmt()
mesh3D.checkConsistencyLight()
- renum = mc.DataArrayInt(60); renum[:15]=range(15,30) ; renum[15:30]=range(15) ; renum[30:45]=range(45,60) ; renum[45:]=range(30,45)
+ renum = mc.DataArrayInt(60)
+ renum[:15]=list(range(15,30))
+ renum[15:30]=list(range(15))
+ renum[30:45]=list(range(45,60))
+ renum[45:]=list(range(30,45))
mesh3D.renumberNodes(renum,60)
Convertir les unités
On vérifie alors que les 3 solutions sont les mêmes : ::
- print cellIdsSol1.getValues()
- print cellIdsSol2.getValues()
- print cellIdsSol3.getValues()
+ print(cellIdsSol1.getValues())
+ print(cellIdsSol2.getValues())
+ print(cellIdsSol3.getValues())
Extraire une sous partie d'un maillage 3D
alors important de voir si ``mesh3DPart`` est bien ordonné, c'est-à-dire si ses cellules sont bien rangées par type géométrique.
On commence par inspecter l'état actuel : ::
- print mesh3DPart.advancedRepr()
+ print(mesh3DPart.advancedRepr())
La fonction suivante fait le même travail : ::
- print mesh3DPart.checkConsecutiveCellTypesAndOrder([mc.NORM_HEXA8, mc.NORM_POLYHED])
+ print(mesh3DPart.checkConsecutiveCellTypesAndOrder([mc.NORM_HEXA8, mc.NORM_POLYHED]))
Ou bien : ::
- print mesh3DPart.checkConsecutiveCellTypes()
+ print(mesh3DPart.checkConsecutiveCellTypes())
On voit que ``mesh3DPart`` contient 6 cellules, quatre HEXA8 puis deux POLYHED. Les cellules sont bien
groupées par type géométrique. Si ce n'était pas le cas, on aurait pu invoquer ``MEDCouplingUMesh.sortCellsInMEDFileFrmt()``.
magn = bary2.magnitude()
ids = magn.findIdsInRange(0.,1e-12)
idStart = int(ids) # ids is assumed to contain only one value, if not an exception is thrown
- ze_range = range(idStart,mesh3D.getNumberOfCells(),mesh2D.getNumberOfCells())
+ ze_range = list(range(idStart,mesh3D.getNumberOfCells(),mesh2D.getNumberOfCells()))
cellIds2Sol2 = extMesh.getMesh3DIds()[ze_range]
Maintenant on construit cette sous partie de ``mesh3D`` en utilisant ``cellIds2Sol1`` ou ``cellIds2Sol2``: ::
import CORBA
orb=CORBA.ORB_init()
ior=orb.object_to_string(ref_m)
- print ior
+ print(ior)
A simple copy/paste in the ParaViS GUI allows to create the source and to have our
mesh rendered on screen.
Début de l'implémentation
~~~~~~~~~~~~~~~~~~~~~~~~~
-Pour commencer l'exercice importer le module ``MEDCoupling``
+Pour commencer l'exercice importer le module ``medcoupling``
et la classe ``MEDCouplingUMeshServant`` du module Python ``MEDCouplingCorba``. ::
- import MEDCoupling as mc
+ import medcoupling as mc
from MEDCouplingCorba import MEDCouplingUMeshServant
Créer un maillage
import CORBA
orb = CORBA.ORB_init()
ior = orb.object_to_string(ref_m)
- print ior
+ print(ior)
Puis, via un copier/coller dans l'IHM ParaViS (Menu "Source -> Para MED Corba Plugin Source"), passer l'IOR.
On voit s'afficher notre maillage.
data=MEDFileData("agitateur.med")
ts=data.getFields()[0].getTimeSteps()
- print ts
+ print(ts)
Get the agitator's mesh (in green) at the time-step (2,-1) (see ts).
To this end use the cell field "DISTANCE_INTERFACE_ELEM_BODY_ELEM_DOM" and select
Sum "torqueOnSkin" using DataArrayDouble.accumulate(). ::
zeTorque=torquePerCellOnSkin.accumulate()
- print "couple = %r N.m"%(zeTorque[2])
+ print("couple = %r N.m"%(zeTorque[2]))
Check the previously computed torque by dividing the power by the angular speed.
Compute the power per skin cell and sum it. ::
speedOnSkin=speedMc.getArray()[tupleIdsInField]
powerSkin=DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
power=powerSkin.accumulate()[0]
- print "power = %r W"%(power)
+ print("power = %r W"%(power))
Compute the angular speed: compute the sum of x^2, y^2 and xz of "posSkin" and build
with NumPy the 2x2 matrix
inertiaSkinValues,inertiaSkinVects=linalg.eig(inertiaSkin)
pos=max(enumerate(inertiaSkinValues),key=lambda x: x[1])[0]
vect0=inertiaSkinVects[pos].tolist()[0]
- print vect0
+ print(vect0)
Thanks to the previous computation we can see that the agitator had a rotation of
1.1183827931 radian (see solution).
Compute and compare the torque on the agitator. ::
omega=1.1183827931/(ts[-1][2]-ts[0][2])
- print "At time-step (%d,%d) at %r s the torque is: %r N.m, power/omega=%r N.m"%(ts[2][0],ts[2][1],ts[2][2],zeTorque[2],power/omega)
+ print("At time-step (%d,%d) at %r s the torque is: %r N.m, power/omega=%r N.m"%(ts[2][0],ts[2][1],ts[2][2],zeTorque[2],power/omega))
Solution
~~~~~~~~
data = ml.MEDFileData("agitateur.med")
ts = data.getFields()[0].getTimeSteps()
- print ts
+ print(ts)
Récupérer le maillage de l'agitateur (en vert) au pas de temps (2,-1) (cf. ts).
La position de l'agitateur est définie par un champ sur le maillage global du système et n'a pas de maillage propre.
Sommer ``torqueOnSkin`` en utilisant la méthode ``DataArrayDouble.accumulate()``. ::
zeTorque = torquePerCellOnSkin.accumulate()
- print "couple = %r N.m" % zeTorque[2]
+ print("couple = %r N.m" % zeTorque[2])
Vérifions le couple calculé précédemment en divisant la puissance par la vitesse *angulaire*.
La vitesse *linéaire* est stockée dans le champ "VITESSE_ELEM_DOM".
speedOnSkin = speedMc.getArray()[tupleIdsInField]
powerSkin = ml.DataArrayDouble.Dot(forceVectSkin,speedOnSkin)
power = powerSkin.accumulate()[0]
- print "power = %r W"%(power)
+ print("power = %r W"%(power))
Calculer la vitesse *angulaire*. Pour ce faire, calculer la somme de ``x^2``, ``y^2`` et ``xz`` de ``posSkin`` et
construire (avec NumPy) la matrice 2x2 d'inertie ``inertiaSkin=[[x2,xy], [xy,z2]]``.
inertiaSkinValues, inertiaSkinVects = np.linalg.eig(inertiaSkin)
pos = max(enumerate(inertiaSkinValues), key=lambda x: x[1])[0]
vect0 = inertiaSkinVects[pos].tolist()[0]
- print vect0
+ print(vect0)
Grâce au calcul précédent on peut déduire que l'agitateur a tourné de 1.1183827931 radian (cf. solution complète pour le
détail - on remet les étapes précédentes dans une fonction que l'on applique sur plusieurs pas de temps).
Calculer et comparer le couple sur l'agitateur. ::
omega = 1.1183827931 / (ts[-1][2]-ts[0][2])
- print "At timestep (%d,%d) (physical time=%r s) the torque is: %r N.m, power/omega=%r N.m " % (ts[2][0],ts[2][1],ts[2][2],zeTorque[2],power/omega)
+ print("At timestep (%d,%d) (physical time=%r s) the torque is: %r N.m, power/omega=%r N.m " % (ts[2][0],ts[2][1],ts[2][2],zeTorque[2],power/omega))
Solution
~~~~~~~~
As we are in nodal connectivity mode it means that common nodes have to merged. This is not the case here.
Merge the nodes closer than 1e-10 and assess the impact on the node count of "fixm". ::
- print "nb of nodes in file : %i"%(fixm.getNumberOfNodes())
+ print("nb of nodes in file : %i"%(fixm.getNumberOfNodes()))
fixm.mergeNodes(1e-10)
- print "nb of non duplicated nodes : %i"%(fixm.getNumberOfNodes())
+ print("nb of non duplicated nodes : %i"%(fixm.getNumberOfNodes()))
Same thing for "Mobile.med" (called "mobm"). Repair it by deleting duplicated nodes. ::
Name this new instance "zone1Mobm", remove all orphan nodes and display. ::
zonesInMobm=mobm.partitionBySpreadZone()
- print "number of zones in mobm : %i"%(len(zonesInMobm))
+ print("number of zones in mobm : %i"%(len(zonesInMobm)))
zone1Mobm=mobm[zonesInMobm[0]]
zone1Mobm.zipCoords()
displayVTK(zone1Mobm,"zone1Mobm.vtu")
all oriented consistently.
To check this let's inspect the areas of the 38 cells of partFixm (variable name "areaPartFixm"). ::
- areaPartFixm=partFixm.getMeasureField(ON_CELLS).getArray()
- print areaPartFixm.getValues()
+ areaPartFixm=partFixm.getMeasureField(isAbs=False).getArray()
+ print(areaPartFixm.getValues())
All values are negative: this MED file doesn't respect the MED file convention.
"partFixm" being mis-oriented and the method MEDCouplingUMesh.Intersect2DMeshes() conserving the orientation, "partFixMob" is also mis-oriented.
To cut long story short, we perform comparison on absolute arrays.
-Check then that the first test check#0 is successful
+Check then that the first test check#0 is successful ::
- areaPartFixm=partFixm.getMeasureField(ON_CELLS).getArray()
+ areaPartFixm=partFixm.getMeasureField(isAbs=False).getArray()
areaPartFixm.abs()
- areaPartFixMob=partFixMob.getMeasureField(ON_CELLS).getArray()
+ areaPartFixMob=partFixMob.getMeasureField(isAbs=False).getArray()
areaPartFixMob.abs()
val1=areaPartFixm.accumulate()[0]
val2=areaPartFixMob.accumulate()[0]
- print "Check #0 %lf == %lf a 1e-8 ? %s"%(val1,val2,str(abs(val1-val2)<1e-8))
+ print("Check #0 %lf == %lf a 1e-8 ? %s"%(val1,val2,str(abs(val1-val2)<1e-8)))
Now check#1. Same spirit as in check#0. ::
- areaZone1Mobm=zone1Mobm.getMeasureField(ON_CELLS).getArray()
+ areaZone1Mobm=zone1Mobm.getMeasureField(isAbs=False).getArray()
areaZone1Mobm.abs()
val3=areaZone1Mobm.accumulate()[0]
ids4=iMob.findIdsNotEqual(-1)
areaPartFixMob2=areaPartFixMob[ids4]
val4=areaPartFixMob2.accumulate()[0]
- print "Check #1 %lf == %lf a 1e-8 ? %s"%(val3,val4,str(abs(val3-val4)<1e-8))
+ print("Check #1 %lf == %lf a 1e-8 ? %s"%(val3,val4,str(abs(val3-val4)<1e-8)))
Finally check#2. ::
isCheck2OK=True
- for icell in xrange(partFixm.getNumberOfCells()):
+ for icell in list(range(partFixm.getNumberOfCells())):
ids5=iPart.findIdsEqual(icell)
areaOfCells=areaPartFixMob[ids5]
areaOfCells.abs()
isCheck2OK=False
pass
pass
- print "Check #2? %s"%(str(isCheck2OK))
+ print("Check #2? %s"%(str(isCheck2OK)))
Use intersection information to create fields
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
en connectivité nodale, il faut absolument que les noeuds soient les mêmes. Il s'avère que cela n'est pas le cas ici.
Fusionner le noeuds distants de moins de 1e-10 et regarder l'impact sur le nombre de noeuds de ``fixm``. ::
- print "Nb of nodes in the file : %i " % (fixm.getNumberOfNodes())
+ print("Nb of nodes in the file : %i " % (fixm.getNumberOfNodes()))
fixm.mergeNodes(1e-10)
- print "Nb of non duplicated nodes : %i" % (fixm.getNumberOfNodes())
+ print("Nb of non duplicated nodes : %i" % (fixm.getNumberOfNodes()))
Même traitement pour ``Mobile.med``, le lire avec l'API avancée de MEDLoader (appeler ``mobm`` l'instance du maillage)
et le réparer en supprimant les noeuds dupliqués. ::
Enfin l'afficher : ::
zonesInMobm = mobm.partitionBySpreadZone()
- print "Nb of zones in mobm : %i" % (len(zonesInMobm))
+ print("Nb of zones in mobm : %i" % (len(zonesInMobm)))
zone1Mobm = mobm[zonesInMobm[0]]
zone1Mobm.zipCoords()
displayVTK(zone1Mobm, "zone1Mobm.vtu")
sont toutes bien orientées ou à minima toutes orientées de la même manière.
Pour ce faire, regardons les aires des 38 cellules de ``partFixm`` (nom de variable : ``areaPartFixm``). ::
- areaPartFixm = partFixm.getMeasureField(ml.ON_CELLS).getArray()
- print areaPartFixm.getValues()
+ areaPartFixm = partFixm.getMeasureField(True).getArray()
+ print(areaPartFixm.getValues())
On voit que toutes les valeurs sont négatives. *Bilan*: ce fichier MED ne respecte pas la convention MED fichier !
``partFixm`` étant mal orienté, et ``MEDCouplingUMesh.Intersect2DMeshes()`` conservant l'orientation,
``partFixMob`` est lui aussi mal orienté.
Bref, on va faire les comparaisons sur des tableaux de valeurs absolues. Vérifier alors **Check #0**. ::
- areaPartFixm = partFixm.getMeasureField(ml.ON_CELLS).getArray()
+ areaPartFixm = partFixm.getMeasureField(isAbs=False).getArray() # prise en compte de l'orientation des mailles
areaPartFixm.abs()
- areaPartFixMob = partFixMob.getMeasureField(ml.ON_CELLS).getArray()
+ areaPartFixMob = partFixMob.getMeasureField(isAbs=False).getArray()
areaPartFixMob.abs()
val1=areaPartFixm.accumulate()[0]
val2=areaPartFixMob.accumulate()[0]
- print "Check #0 %lf == %lf with precision 1e-8? %s" % (val1,val2,str(abs(val1-val2)<1e-8))
+ print("Check #0 %lf == %lf with precision 1e-8? %s" % (val1,val2,str(abs(val1-val2)<1e-8)))
On peut passer au **Check #1**. L'esprit est le même que le **Check #0**. ::
- areaZone1Mobm = zone1Mobm.getMeasureField(ml.ON_CELLS).getArray()
+ areaZone1Mobm = zone1Mobm.getMeasureField(isAbs=False).getArray()
areaZone1Mobm.abs()
val3 = areaZone1Mobm.accumulate()[0]
ids4 = iMob.findIdsNotEqual(-1)
areaPartFixMob2 = areaPartFixMob[ids4]
val4 = areaPartFixMob2.accumulate()[0]
- print "Check #1 %lf == %lf with precision 1e-8 ? %s" % (val3,val4,str(abs(val3-val4)<1e-8))
+ print("Check #1 %lf == %lf with precision 1e-8 ? %s" % (val3,val4,str(abs(val3-val4)<1e-8)))
Puis le **Check #2**. ::
isCheck2OK = True
- for icell in xrange(partFixm.getNumberOfCells()):
+ for icell in list(range(partFixm.getNumberOfCells())):
ids5 = iPart.findIdsEqual(icell)
areaOfCells = areaPartFixMob[ids5]
areaOfCells.abs()
isCheck2OK = False
pass
pass
- print "Check #2? %s" % (str(isCheck2OK))
+ print("Check #2? %s" % (str(isCheck2OK)))
Utiliser les informations de l'intersection pour en faire des champs
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Début de l'implémentation
~~~~~~~~~~~~~~~~~~~~~~~~~
-Pour commencer l'exercice importer le module Python ``MEDCoupling``, ``MEDCouplingRemapper``, ``numpy``, ``scipy``, ``multiprocessing``
+Pour commencer l'exercice importer le module Python ``medcoupling``, ``MEDCouplingRemapper``, ``numpy``, ``scipy``, ``multiprocessing``
et ``datetime`` pour chronométrer : ::
- import MEDCoupling as mc
+ import medcoupling as mc
import MEDCouplingRemapper as mr
import multiprocessing as mp
from datetime import datetime
remap=mr.MEDCouplingRemapper()
strt=datetime.now()
assert(remap.prepare(m,m2,"P0P0")==1)
- print "time in sequential : %s"%(str(datetime.now()-strt))
+ print("time in sequential : %s"%(str(datetime.now()-strt)))
Stockons la sparse matrix scipy dans ``matSeq``. ::
On peut se faire aider de ``mc.DataArray.GetSlice``. ::
workToDo=[]
- for i in xrange(nbProc):
+ for i in list(range(nbProc)):
s=mc.DataArray.GetSlice(slice(0,m2.getNumberOfCells(),1),i,nbProc)
part=m2[s]
partToGlob=mc.DataArrayInt.Range(s.start,s.stop,s.step)
pool = mp.Pool()
asyncResult = pool.map_async(work,workToDo)
subMatrices = asyncResult.get()
- print "time in parallel (x%d) : %s"%(nbProc,str(datetime.now()-strt))
+ print("time in parallel (x%d) : %s"%(nbProc,str(datetime.now()-strt)))
.. note:: A noter la magie ! On a transféré entre le process maitre et les process esclave sans même s'en rendre compte les maillages et les DataArrayInt contenus dans ``workToDo`` !
Merci à la pickelisation des objets MEDCoupling :)
Début de l'implémentation
~~~~~~~~~~~~~~~~~~~~~~~~~
-Pour commencer l'exercice importer le module Python ``MEDCoupling``: ::
+Pour commencer l'exercice importer le module Python ``medcoupling``: ::
- import MEDCoupling as mc
+ import medcoupling as mc
NumPy est un prérequis optionnel, vérifions que nous en bénéficions bien : ::
Et afficher ``nparr``. ::
- print nparr.__repr__()
- print nparr.tolist()
+ print(nparr.__repr__())
+ print(nparr.tolist())
Mais est ce qu'on ne nous a pas mystifié ? ``arr`` et ``nparr`` partagent-ils le même bloc mémoire ?
Pour le vérifier assignons 7.0 un tuple sur 2 avec ``nparr`` et vérifions que ``arr`` et ``nparr`` sont simultanément modifiés. ::
nparr[::2] = 7.
- print nparr.__repr__()
- print arr.__repr__()
+ print(nparr.__repr__())
+ print(arr.__repr__())
C'est rigolo ! Mais si je détruis ``arr`` (le premier à avoir alloué la mémoire) est-ce que ``nparr`` est tué aussi ?
Ne risque-t-on pas le SIGSEGV ?
del arr
import gc; gc.collect() # Make sure the object has been deleted
- print nparr.__repr__()
+ print(nparr.__repr__())
OK super. Mais inversement puis je faire une instance de ``DataArrayDouble`` avec ``nparr`` ? Oui, en utilisant le constructeur
qui prend un ``nparray`` en entrée.
Et afficher le contenu.::
arr2 = mc.DataArrayDouble(nparr)
- print arr2.__repr__()
+ print(arr2.__repr__())
Modifions ``nparr`` en assignant 5.0 pour tous les tuples et vérifier que les 2 représentations ont bien été modifiées simultanément.::
nparr[:] = 5.
- print nparr.__repr__()
- print arr2.__repr__()
+ print(nparr.__repr__())
+ print(arr2.__repr__())
Nous en profitons pour montrer un petit service pratique avec NumPy, à savoir, l'écriture optimisée.
Ecrivons le contenu binaire de ``nparr`` dans un fichier. ::
a[:] = 3.14
f.flush()
b = np.memmap(f2,dtype='float64',mode='r',offset=0,shape=(12,))
- print b.__repr__()
+ print(b.__repr__())
On voit donc que le passage de MEDCoupling à NumPy se fait directement et de manière optimisée. Donc ca peut valoir le coup !
Tout ce qui vient d'être montré marche aussi avec des ``DataArrayInt``.
Début de l'implémentation
~~~~~~~~~~~~~~~~~~~~~~~~~
-Pour commencer l'exercice importer le module Python ``MEDCoupling``: ::
+Pour commencer l'exercice importer le module Python ``medcoupling``: ::
- import MEDCoupling as mc
+ import medcoupling as mc
NumPy est un prérequis optionnel, vérifions que nous en bénéficions bien : ::
Et afficher ``nparr``. ::
- print nparr.__repr__()
- print nparr.tolist()
+ print(nparr.__repr__())
+ print(nparr.tolist())
Mais est ce qu'on ne nous a pas mystifié ? ``arr`` et ``nparr`` partagent-ils le même bloc mémoire ?
Pour le vérifier assignons 7.0 un tuple sur 2 avec ``nparr`` et vérifions que ``arr`` et ``nparr`` sont simultanément modifiés. ::
nparr[::2] = 7.
- print nparr.__repr__()
- print arr.__repr__()
+ print(nparr.__repr__())
+ print(arr.__repr__())
C'est rigolo ! Mais si je détruis ``arr`` (le premier à avoir alloué la mémoire) est-ce que ``nparr`` est tué aussi ?
Ne risque-t-on pas le SIGSEGV ?
del arr
import gc; gc.collect() # Make sure the object has been deleted
- print nparr.__repr__()
+ print(nparr.__repr__())
OK super. Mais inversement puis je faire une instance de ``DataArrayDouble`` avec ``nparr`` ? Oui, en utilisant le constructeur
qui prend un ``nparray`` en entrée.
Et afficher le contenu.::
arr2 = mc.DataArrayDouble(nparr)
- print arr2.__repr__()
+ print(arr2.__repr__())
Modifions ``nparr`` en assignant 5.0 pour tous les tuples et vérifier que les 2 représentations ont bien été modifiées simultanément.::
nparr[:] = 5.
- print nparr.__repr__()
- print arr2.__repr__()
+ print(nparr.__repr__())
+ print(arr2.__repr__())
Nous en profitons pour montrer un petit service pratique avec NumPy, à savoir, l'écriture optimisée.
Ecrivons le contenu binaire de ``nparr`` dans un fichier. ::
a[:] = 3.14
f.flush()
b = np.memmap(f2,dtype='float64',mode='r',offset=0,shape=(12,))
- print b.__repr__()
+ print(b.__repr__())
On voit donc que le passage de MEDCoupling à NumPy se fait directement et de manière optimisée. Donc ca peut valoir le coup !
Tout ce qui vient d'être montré marche aussi avec des ``DataArrayInt``.
::
myMatrix=remap.getCrudeMatrix()
- print myMatrix # to see what it looks like
+ print(myMatrix) # to see what it looks like
sumByRows=DataArrayDouble(len(myMatrix))
for i,wIt in enumerate(sumByRows):
su=0.
for it in myMatrix[i]:
su+=myMatrix[i][it]
wIt[0]=su
- print "Does interpolation look OK? %s"%(str(sumByRows.isUniform(1.,1e-12)))
+ print("Does interpolation look OK? %s"%(str(sumByRows.isUniform(1.,1e-12))))
.. note:: Some triangles were added into "srcMesh" to make "myMatrix" less boring. "myMatrix".
Check that with this nature the field integral is conserved. On the other side
the sum on cells (accumulation) is NOT conserved. ::
- print "IntensiveMaximum %lf == %lf"%(srcField.integral(True)[0],trgFieldCV.integral(True)[0])
- print "IntensiveMaximum %lf != %lf"%(srcField.getArray().accumulate()[0],trgFieldCV.getArray().accumulate()[0])
+ print("IntensiveMaximum %lf == %lf"%(srcField.integral(True)[0],trgFieldCV.integral(True)[0]))
+ print("IntensiveMaximum %lf != %lf"%(srcField.getArray().accumulate()[0],trgFieldCV.getArray().accumulate()[0]))
Set the nature of "srcField" to ExtensiveConservation (extensive field, e.g. a power). ::
cumulative sum on cells is conserved. ::
::
- print "ExtensiveConservation %lf != %lf"%(srcField.integral(True)[0],trgFieldI.integral(True)[0])
- print "ExtensiveConservation %lf == %lf"%(srcField.getArray().accumulate()[0],trgFieldI.getArray().accumulate()[0])
+ print("ExtensiveConservation %lf != %lf"%(srcField.integral(True)[0],trgFieldI.integral(True)[0]))
+ print("ExtensiveConservation %lf == %lf"%(srcField.getArray().accumulate()[0],trgFieldI.getArray().accumulate()[0]))
Visualize the fields using ParaViS.
Début de l'implémentation
~~~~~~~~~~~~~~~~~~~~~~~~~
-Pour commencer l'exercice importer le module Python ``MEDCoupling`` et la
+Pour commencer l'exercice importer le module Python ``medcoupling`` et la
classe ``MEDCouplingRemapper`` du module ``MEDCouplingRemapper``. ::
- import MEDCoupling as mc
+ import medcoupling as mc
from MEDCouplingRemapper import MEDCouplingRemapper
Vérifier notamment que pour chaque cellule de ``trgMesh`` la somme des aires fait toujours 1. ::
myMatrix = remap.getCrudeMatrix()
- print myMatrix
+ print(myMatrix)
sumByRows = mc.DataArrayDouble(len(myMatrix))
for i,wIt in enumerate(sumByRows):
su = 0.
for it in myMatrix[i]:
su += myMatrix[i][it]
wIt[0] = su
- print "Is interpolation well prepared?", sumByRows.isUniform(1.,1e-12)
+ print("Is interpolation well prepared?", sumByRows.isUniform(1.,1e-12))
.. note:: Les triangles dans ``srcMesh`` ont été rajoutés pour casser la monotonie de la matrice ``myMatrix``.
integSource = srcField.integral(True)[0]
integTarget = trgFieldCV.integral(True)[0]
- print "IntensiveMaximum -- integrals: %lf == %lf" % (integSource, integTarget)
+ print("IntensiveMaximum -- integrals: %lf == %lf" % (integSource, integTarget))
accSource = srcField.getArray().accumulate()[0]
accTarget = trgFieldCV.getArray().accumulate()[0]
- print "IntensiveMaximum -- sums: %lf != %lf" % (accSource, accTarget)
+ print("IntensiveMaximum -- sums: %lf != %lf" % (accSource, accTarget))
Maintenant mettre la nature de ``srcField`` à ``ExtensiveConservation``. Le champ doit être interprété commé étant
integSource = srcField.integral(True)[0]
integTarget = trgFieldI.integral(True)[0]
- print "ExtensiveConservation -- integrals: %lf != %lf" % (integSource, integTarget)
+ print("ExtensiveConservation -- integrals: %lf != %lf" % (integSource, integTarget))
accSource = srcField.getArray().accumulate()[0]
accTarget = trgFieldI.getArray().accumulate()[0]
- print "ExtensiveConservation -- sums: %lf == %lf" % (accSource, accTarget)
+ print("ExtensiveConservation -- sums: %lf == %lf" % (accSource, accTarget))
Visualiser les champs avec ParaViS, ou en les écrivant dans un fichier.
NodeField1=NodeField[proc1] ; CellField1=CellField[proc1] ; CellField1.setMesh(NodeField1.getMesh())
proc0_fname="proc0.med"
- MEDLoader.WriteField(proc0_fname,NodeField0,True)
- MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc0_fname,CellField0)
+ WriteField(proc0_fname,NodeField0,True)
+ WriteFieldUsingAlreadyWrittenMesh(proc0_fname,CellField0)
proc1_fname="proc1.med"
- MEDLoader.WriteField(proc1_fname,NodeField1,True)
- MEDLoader.WriteFieldUsingAlreadyWrittenMesh(proc1_fname,CellField1)
+ WriteField(proc1_fname,NodeField1,True)
+ WriteFieldUsingAlreadyWrittenMesh(proc1_fname,CellField1)
Reading and merging 2 MED files - Easy (but non optimal) version
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In the two files "proc0.med" and "proc1.med" read the respective "CellField" with the basic API. Aggregate the two and store the result in "CellField_read". ::
- CellField0_read=MEDLoader.ReadFieldCell("proc0.med","mesh",0,"CellField",5,6)
- CellField1_read=MEDLoader.ReadFieldCell("proc1.med","mesh",0,"CellField",5,6)
+ CellField0_read=ReadFieldCell("proc0.med","mesh",0,"CellField",5,6)
+ CellField1_read=ReadFieldCell("proc1.med","mesh",0,"CellField",5,6)
CellField_read=MEDCouplingFieldDouble.MergeFields([CellField0_read,CellField1_read])
.. note:: It might seem to the reader that the cell type information is repeated uselessly, but don't forget that the MED file norm doesn't forbid a field to be defined simultaneously on nodes and on Gauss points for example ...
CellFieldCpy=CellField.deepCopy()
CellFieldCpy.substractInPlaceDM(CellField_read,10,1e-12)
CellFieldCpy.getArray().abs()
- print CellFieldCpy.getArray().isUniform(0.,1e-12)
+ print(CellFieldCpy.getArray().isUniform(0.,1e-12))
Let's do the same on "NodeField". The main difference here is that redundant information is created at the boundary. ::
- NodeField0_read=MEDLoader.ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
- NodeField1_read=MEDLoader.ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
+ NodeField0_read=ReadFieldNode("proc0.med","mesh",0,"NodeField",5,6)
+ NodeField1_read=ReadFieldNode("proc1.med","mesh",0,"NodeField",5,6)
NodeField_read=MEDCouplingFieldDouble.MergeFields([NodeField0_read,NodeField1_read])
.. note:: The mesh is read a second time here, which can be damaging in terms of performance.
Compare "NodeFieldCpy" and "NodeField_read" still using MEDCouplingFieldDouble.substractInPlaceDM(). ::
NodeFieldCpy.substractInPlaceDM(NodeField_read,10,1e-12)
- print NodeFieldCpy.getArray().isUniform(0.,1e-12)
+ print(NodeFieldCpy.getArray().isUniform(0.,1e-12))
Read/write of two separated MED files - More complex but more efficient version
cellFieldCpy = cellField.deepCopy()
cellFieldCpy.substractInPlaceDM(cellField_read,10,1e-12)
cellFieldCpy.getArray().abs()
- print cellFieldCpy.getArray().isUniform(0.,1e-12)
+ print(cellFieldCpy.getArray().isUniform(0.,1e-12))
Opérons le même travail sur "NodeField" que celui réalisé plus haut sur "CellField".
La différence ici c'est qu'il va y avoir duplication de l'information à la frontière, car les noeuds limites sont partagés
Comparer ``nodeFieldCpy`` et ``nodeField_read`` toujours en utilisant ``MEDCouplingFieldDouble.substractInPlaceDM()`` : ::
nodeFieldCpy.substractInPlaceDM(nodeField_read,10,1e-12)
- print nodeFieldCpy.getArray().isUniform(0.,1e-12)
+ print(nodeFieldCpy.getArray().isUniform(0.,1e-12))
Lecture et merge des 2 fichiers MED séparés (moins facile, mais plus optimal)
~~~~~~~~~~~~~~~~~~~~
To implement this exercise we use the Python scripting language and import the MEDLoader Python module.
-The whole MEDCoupling module is fully included in MEDLoader. No need to import MEDCoupling when MEDLoader has been loaded. ::
+The whole MEDCoupling module is fully included in MEDLoader. No need to import medcoupling when MEDLoader has been loaded. ::
import MEDLoader as ml
targetConn=[0,3,4,1, 1,4,2, 4,5,2, 6,7,4,3, 7,8,5,4]
targetMesh=ml.MEDCouplingUMesh.New("MyMesh",2)
targetMesh.allocateCells(5)
- targetMesh.insertNextCell(NORM_TRI3,3,targetConn[4:7])
- targetMesh.insertNextCell(NORM_TRI3,3,targetConn[7:10])
- targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[0:4])
- targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[10:14])
- targetMesh.insertNextCell(NORM_QUAD4,4,targetConn[14:18])
+ targetMesh.insertNextCell(ml.NORM_TRI3,3,targetConn[4:7])
+ targetMesh.insertNextCell(ml.NORM_TRI3,3,targetConn[7:10])
+ targetMesh.insertNextCell(ml.NORM_QUAD4,4,targetConn[0:4])
+ targetMesh.insertNextCell(ml.NORM_QUAD4,4,targetConn[10:14])
+ targetMesh.insertNextCell(ml.NORM_QUAD4,4,targetConn[14:18])
myCoords=ml.DataArrayDouble.New(targetCoords,9,2)
targetMesh.setCoords(myCoords)
Then we are ready to write targetMesh and targetMesh1 into TargetMesh2.med. ::
- meshMEDFile=MEDFileUMesh.New()
+ meshMEDFile=ml.MEDFileUMesh.New()
meshMEDFile.setMeshAtLevel(0,targetMesh)
meshMEDFile.setMeshAtLevel(-1,targetMesh1)
meshMEDFile.write("TargetMesh2.med",2) # 2 stands for write from scratch
-Create 2 groups on level 0. The first called "grp0_Lev0" on cells [0,1,3] and the second called "grp1_Lev0" on cells [1,2,3,4] ::
+Create 2 groups on level 0. The first called "grp0_Lev0" on cells [0,1,3] and the second called "grp1_Lev0" on cells [1,2,3,4] ::
grp0_0=ml.DataArrayInt.New([0,1,3]) ; grp0_0.setName("grp0_Lev0")
grp1_0=ml.DataArrayInt.New([1,2,3,4]) ; grp1_0.setName("grp1_Lev0")
Then trying to read it. ::
- meshMEDFileRead=MEDFileMesh.New("TargetMesh2.med")
+ meshMEDFileRead=ml.MEDFileMesh.New("TargetMesh2.med")
meshRead0=meshMEDFileRead.getMeshAtLevel(0)
meshRead1=meshMEDFileRead.getMeshAtLevel(-1)
- print "Is the mesh at level 0 read in file equals targetMesh ? %s"%(meshRead0.isEqual(targetMesh,1e-12))
- print "Is the mesh at level -1 read in file equals targetMesh ? %s"%(meshRead1.isEqual(targetMesh1,1e-12))
+ print("Is the mesh at level 0 read in file equals targetMesh ? %s"%(meshRead0.isEqual(targetMesh,1e-12)))
+ print("Is the mesh at level -1 read in file equals targetMesh ? %s"%(meshRead1.isEqual(targetMesh1,1e-12)))
Print available levels for group "grp0_Lev0" ::
- print meshMEDFileRead.getGrpNonEmptyLevels("grp0_Lev0")
+ print(meshMEDFileRead.getGrpNonEmptyLevels("grp0_Lev0"))
Request for cell ids of group "grp0_Lev0" ::
grp0_0_read=meshMEDFileRead.getGroupArr(0,"grp0_Lev0")
- print "Is group \"grp0_Lev0\" are the same ? %s"%(grp0_0_read.isEqual(grp0_0))
+ print("Is group \"grp0_Lev0\" are the same ? %s"%(grp0_0_read.isEqual(grp0_0)))
Writing and Reading fields
~~~~~~~~~~~~~~~~~~~~~~~~~~
Creation of a simple vector field on cells called f. ::
- f=ml.MEDCouplingFieldDouble.New(ON_CELLS,ONE_TIME)
+ f=ml.MEDCouplingFieldDouble.New(ml.ON_CELLS,ml.ONE_TIME)
f.setTime(5.6,7,8)
f.setArray(targetMesh.computeCellCenterOfMass())
f.setMesh(targetMesh)
Put f into a MEDFileField1TS for preparation of MED writing ::
- fMEDFile=MEDFileField1TS.New()
+ fMEDFile=ml.MEDFileField1TS.New()
fMEDFile.setFieldNoProfileSBT(f)
Append field to "TargetMesh2.med" ::
Read it : ::
- fMEDFileRead=MEDFileField1TS.New("TargetMesh2.med",f.getName(),7,8)
- fRead1=fMEDFileRead.getFieldOnMeshAtLevel(ON_CELLS,0,meshMEDFileRead) # fastest method. No reading of the supporting mesh.
- fRead2=fMEDFileRead.getFieldAtLevel(ON_CELLS,0) # like above but mesh is re-read from file...
- print "Does the field f remain the same using fast method ? %s"%(fRead1.isEqual(f,1e-12,1e-12))
- print "Does the field f remain the same using slow method ? %s"%(fRead2.isEqual(f,1e-12,1e-12))
+ fMEDFileRead=ml.MEDFileField1TS.New("TargetMesh2.med",f.getName(),7,8)
+ fRead1=fMEDFileRead.getFieldOnMeshAtLevel(ml.ON_CELLS,0,meshMEDFileRead) # fastest method. No reading of the supporting mesh.
+ fRead2=fMEDFileRead.getFieldAtLevel(ml.ON_CELLS,0) # like above but mesh is re-read from file...
+ print("Does the field f remain the same using fast method ? %s"%(fRead1.isEqual(f,1e-12,1e-12)))
+ print("Does the field f remain the same using slow method ? %s"%(fRead2.isEqual(f,1e-12,1e-12)))
Writing and Reading fields on a "profile"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Put it into MEDFileField1TS data structure. ::
- fMEDFile2=MEDFileField1TS.New()
+ fMEDFile2=ml.MEDFileField1TS.New()
fMEDFile2.setFieldProfile(fPart,meshMEDFileRead,0,pfl)
fMEDFile2.write("TargetMesh2.med",0) # 0 is very important here because we want to append to TargetMesh2.med and not to scratch it
Read "fPart" field from File "TargetMesh2.med". ::
- fMEDFileRead2=MEDFileField1TS.New("TargetMesh2.med",fPart.getName(),7,8)
- fPartRead,pflRead=fMEDFileRead2.getFieldWithProfile(ON_CELLS,0,meshMEDFileRead)
- print fPartRead.isEqualWithoutConsideringStr(fPart.getArray(),1e-12)
- print pflRead.isEqualWithoutConsideringStr(pfl)
+ fMEDFileRead2=ml.MEDFileField1TS.New("TargetMesh2.med",fPart.getName(),7,8)
+ fPartRead,pflRead=fMEDFileRead2.getFieldWithProfile(ml.ON_CELLS,0,meshMEDFileRead)
+ print(fPartRead.isEqualWithoutConsideringStr(fPart.getArray(),1e-12))
+ print(pflRead.isEqualWithoutConsideringStr(pfl))
Solution
~~~~~~~~
meshMEDFileRead = ml.MEDFileMesh.New("TargetMesh2.med") # a new is needed because it returns a MEDFileUMesh (MEDFileMesh is abstract)
meshRead0 = meshMEDFileRead.getMeshAtLevel(0)
meshRead1 = meshMEDFileRead.getMeshAtLevel(-1)
- print "Is level 0 in the file equal to 'targetMesh'?", meshRead0.isEqual(targetMesh,1e-12)
- print "Is level 0 in the file equal to 'targetMesh1'?", meshRead1.isEqual(targetMesh1,1e-12)
+ print("Is level 0 in the file equal to 'targetMesh'?", meshRead0.isEqual(targetMesh,1e-12))
+ print("Is level 0 in the file equal to 'targetMesh1'?", meshRead1.isEqual(targetMesh1,1e-12))
Affichons les niveaux disponibles pour le groupe ``grp0_Lev0`` : ::
- print meshMEDFileRead.getGrpNonEmptyLevels("grp0_Lev0")
+ print(meshMEDFileRead.getGrpNonEmptyLevels("grp0_Lev0"))
Et récupérons enfin les identifiants de cellules contenus dans le groupe ``grp0_Lev0`` : ::
grp0_0_read = meshMEDFileRead.getGroupArr(0,"grp0_Lev0")
- print "Is group 'grp0_Lev0' equal to what is read in the file?" , grp0_0_read.isEqual(grp0_0)
+ print("Is group 'grp0_Lev0' equal to what is read in the file?" , grp0_0_read.isEqual(grp0_0))
Lire/écrire des champs avec l'API avancée
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
fMEDFileRead = ml.MEDFileField1TS("TargetMesh2.med",f.getName(),7,8)
fRead1 = fMEDFileRead.getFieldOnMeshAtLevel(ml.ON_CELLS,0,meshMEDFileRead) # Quickest way, not re-reading mesh in the file.
fRead2 = fMEDFileRead.getFieldAtLevel(ml.ON_CELLS,0) # Like above, but this time the mesh is read!
- print "Does the field remain OK with the quick method?", fRead1.isEqual(f,1e-12,1e-12)
- print "Does the field remain OK with the slow method?", fRead2.isEqual(f,1e-12,1e-12)
+ print("Does the field remain OK with the quick method?", fRead1.isEqual(f,1e-12,1e-12))
+ print("Does the field remain OK with the slow method?", fRead2.isEqual(f,1e-12,1e-12))
Lire/écrire un champ sur un "profil"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
fMEDFileRead2 = ml.MEDFileField1TS("TargetMesh2.med",fPart.getName(),7,8)
fPartRead, pflRead = fMEDFileRead2.getFieldWithProfile(ml.ON_CELLS,0,meshMEDFileRead)
- print "Is the partial field correctly read?", fPartRead.isEqualWithoutConsideringStr(fPart.getArray(),1e-12)
- print "Is the list of cell identifiers matching?", pflRead.isEqualWithoutConsideringStr(pfl)
+ print("Is the partial field correctly read?", fPartRead.isEqualWithoutConsideringStr(fPart.getArray(),1e-12))
+ print("Is the list of cell identifiers matching?", pflRead.isEqualWithoutConsideringStr(pfl))
Solution
~~~~~~~~
Then trying to read it. ::
meshRead=ml.ReadUMeshFromFile("TargetMesh.med",targetMesh.getName(),0)
- print "Is the mesh read in file equals targetMesh? %s"%(meshRead.isEqual(targetMesh,1e-12))
+ print("Is the mesh read in file equals targetMesh? %s"%(meshRead.isEqual(targetMesh,1e-12)))
Writing/Reading a field on one time step at once
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Reading into MyFirstField.med ::
f2=ml.ReadFieldCell("MyFirstField.med",f.getMesh().getName(),0,f.getName(),7,8)
- print "Is the field read in file equals f ? %s"%(f2.isEqual(f,1e-12,1e-12))
+ print("Is the field read in file equals f ? %s"%(f2.isEqual(f,1e-12,1e-12)))
Writing/Reading a field on one or many times steps in "multi-session mode"
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
... et relu. ::
meshRead = ml.ReadUMeshFromFile("TargetMesh.med",targetMesh.getName(),0)
- print "Is the read mesh equal to 'targetMesh' ?", meshRead.isEqual(targetMesh,1e-12)
+ print("Is the read mesh equal to 'targetMesh' ?", meshRead.isEqual(targetMesh,1e-12))
Lire/Ecrire un champ sur un pas de temps
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Nous relisons ensuite MyFirstField.med : ::
f2 = ml.ReadFieldCell("MyFirstField.med", f.getMesh().getName(), 0, f.getName(), 7, 8)
- print "Is the read field identical to 'f' ?", f2.isEqual(f,1e-12,1e-12)
+ print("Is the read field identical to 'f' ?", f2.isEqual(f,1e-12,1e-12))
.. note:: Lors de la lecture du champ, on doit donc connaître: son nom, le nom de sa mesh de support
et le pas de temps voulu. Des fonctions du type ``MEDFileFields.getFieldsNames()`` ou encore
Nous pouvons relire tout cela avec des méthodes similaires à ce qui a été vu précédemment : ::
f3 = ml.ReadFieldCell("MySecondField.med",f.getMesh().getName(),0,f.getName(),7,8)
- print "Is the field read in file equals to 'f' ?", f.isEqual(f3,1e-12,1e-12)
+ print("Is the field read in file equals to 'f' ?", f.isEqual(f3,1e-12,1e-12))
f4 = ml.ReadFieldCell("MySecondField.med",f.getMesh().getName(),0,f.getName(),9,10)
- print "Is the field read in file equals to 'f2' ?", f2.isEqual(f4,1e-12,1e-12)
+ print("Is the field read in file equals to 'f2' ?", f2.isEqual(f4,1e-12,1e-12))
Solution
~~~~~~~~
ipar.append("AP_MODULES_LIST", "Mesh")
import sys
-if sys.platform == "win32":
- from MEDCouplingCompat import *
-else:
- from MEDCoupling import *
+from medcoupling import *
from MEDLoader import WriteMesh
-arr=DataArrayDouble(range(2))
+arr=DataArrayDouble(list(range(2)))
mc=MEDCouplingCMesh("m1")
mc.setCoords(arr,arr)
m1=mc.buildUnstructured()
smesh = smeshBuilder.New(salome.myStudy)
([mesh_1,mesh_2], status) = smesh.CreateMeshesFromMED(r'mesh1.med')
-g1 = mesh_1.MakeGroupByIds("all nodes", SMESH.NODE, range(1,10))
-g2 = mesh_2.MakeGroupByIds("all nodes", SMESH.NODE, range(1,10))
+g1 = mesh_1.MakeGroupByIds("all nodes", SMESH.NODE, list(range(1,10)))
+g2 = mesh_2.MakeGroupByIds("all nodes", SMESH.NODE, list(range(1,10)))
### Store presentation parameters of displayed objects
import iparameters
ipar.append("AP_MODULES_LIST", "Mesh")
import sys
-if sys.platform == "win32":
- from MEDCouplingCompat import *
-else:
- from MEDCoupling import *
+from medcoupling import *
-coordsArr=DataArrayDouble(range(5))
+coordsArr=DataArrayDouble(list(range(5)))
mesh1=MEDCouplingCMesh("mesh")
mesh1.setCoords(coordsArr,coordsArr)
mesh1 = mesh1.buildUnstructured()
medfile2="mesh2.med"
m4=MEDCouplingCMesh("box")
-coo=DataArrayDouble(range(7))
+coo=DataArrayDouble(list(range(7)))
m4.setCoords(coo[:5],coo[:5],coo)
m4=m4.buildUnstructured()
valsArr1=m4.computeCellCenterOfMass()
