Manipuler des champs de double
------------------------------
Les champs dans MEDCoupling ont comme support un unique maillage, de dimension fixĂŠe, et bien dĂŠfini.
Cela semble trivial mais c'est en fait une diffĂŠrence majeure avec la notion de champ dans MED fichier, qui elle est beaucoup
plus permissive.
Les champs sont utiles pour :
* stocker des valeurs d'une grandeur physique relative au problème traitÊ, mais aussi
* des services de haut niveau oĂš l'interaction avec les maillages
est requise comme par exemple ``getValueOn()``, ``getValueOnMulti()``, ``integral()``, ``getMeasureField``
``normL1()``, ``normL2()``, ``fillFromAnalytic()``, ... qui calculent toutes des valeurs en lien avec le maillage
(par exemple le *volume* des cellules)
* expliciter prĂŠcisĂŠment les informations ĂŠchangĂŠes entre les diffĂŠrents codes
lors de couplage.
Pour information, l'implĂŠmentation de ``MEDCouplingFieldDouble`` est relativement petite car cette classe
dÊlègue la très large majoritÊ de ses traitements à des instances de classes aggrÊgÊes
comme ``MEDCouplingMesh``, ``DataArrayDouble``, et ``MEDCouplingSpatialDiscretization``.
La classe ``MEDCouplingFieldDouble`` permet d'assurer la cohĂŠrence entre tous ces ĂŠlĂŠments.
Il est souvent possible et mĂŞme parfois recommandĂŠ de manipuler les tableaux (un ``DataArrayDouble``)
et/ou le maillage d'une instance de ``MEDCouplingFieldDouble`` directement.
Objectifs
~~~~~~~~~
Cet exercice met l'accent sur la relation entre le maillage et les valeurs d'un champ.
* CrĂŠer un champ
* AgrĂŠger des champs
* Construire une sous-partie d'un champ
* RenumĂŠroter les entitĂŠs d'un champ
* Comparer 2 champs venant de 2 sources diffĂŠrentes
* Evaluation d'un champ sur un ensemble de points
* Exploser un champ
DĂŠbut de l'implementation
~~~~~~~~~~~~~~~~~~~~~~~~~
Importer le module Python ``MEDCoupling``. ::
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. ::
xarr = mc.DataArrayDouble.New(11,1)
xarr.iota(0.) # Generate s, s+1, s+2, ... with a given start value s
cmesh = mc.MEDCouplingCMesh.New()
cmesh.setCoords(xarr,xarr,xarr)
mesh = cmesh.buildUnstructured()
.. note:: La mÊthode ``MEDCouplingMesh.buildUnstructured()`` est très utile pour construire rapidement un maillage
non structurĂŠ afin de tester quelque chose.
Afin de mettre en Êvidence le problème des types gÊomÊtriques multiples, convertir en polyhèdres
les cellules d'identifiant pair ::
mesh.convertToPolyTypes(mc.DataArrayInt.Range(0,mesh.getNumberOfCells(),2))
CrĂŠation d'un champ
~~~~~~~~~~~~~~~~~~~
CrĂŠer un champ scalaire (une seule composante) aux cellules (P0) appelĂŠ "MyField" en appliquant la fonction analytique
suivante ``(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)``, oĂš ``(x, y, z)`` reprĂŠsente implicitement les coordonnĂŠes du barycentre
d'une cellule.
Pour cela, deux possiblitĂŠs :
* Directement en appelant ``fillFromAnalytic()`` sur un maillage ::
f = mesh.fillFromAnalytic(mc.ON_CELLS,1,"(x-5.)*(x-5.)+(y-5.)*(y-5.)+(z-5.)*(z-5.)") # 1 means that the field should have one component
f.setName("MyField")
* Ou en crĂŠant au prĂŠalable un champ non initialisĂŠ, et en appliquant ``fillFromAnalytic()`` sur cette
instance de ``MEDCouplingFieldDouble`` ::
f2 = mc.MEDCouplingFieldDouble(mc.ON_CELLS, mc.ONE_TIME)
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
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)
.. note:: Le ``WithoutConsideringStr`` dans le nom de la mĂŠthode prĂŠcĂŠdente indique que les noms des champs ne seront
pas comparĂŠs. On retrouve ce suffixe dans d'autres mĂŠthodes MEDCoupling.
Construire une sous-partie d'un champ
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
RĂŠcupĂŠrer dans une variable ``ids1`` la liste des identifiants de cellules pour lesquelles la valeur du champ est dans le
range [0.0,5.0]. Utiliser pour cela la mĂŠthode ``DataArrayDouble.findIdsInRange()``. Avec ce rĂŠsultat, construire la
sous-partie ``fPart1`` du champ ``f``. ::
da1 = f.getArray() # a DataArrayDouble, which is a direct reference (not a copy) of the field's values
ids1 = da1.findIdsInRange(0., 5.)
fPart1 = f.buildSubPart(ids1)
fPart1.writeVTK("ExoField_fPart1.vtu")
.. image:: images/FieldDouble1_1.png
SĂŠlectionner la partie ``fPart2`` du champ ``f`` dont toutes les valeurs de tuples
sont dans ``[50.,+infinity)``. ::
ids2 = f.getArray().findIdsInRange(50., 1.e300)
fPart2 = f.buildSubPart(ids2)
Ce genre de technique permet d'extraire facilement les parties d'un champ relatives Ă un groupe de mailles par exemple.
RenumĂŠroter les entitĂŠs d'un champ
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
La partie ``fPart1`` gĂŠnĂŠrĂŠe est valide d'un point de vue de MEDCoupling. Mais elle
n'est pas valide d'un point de vue de MED *fichier*.
Une renumÊrotation s'impose dans l'hypothèse de stocker ce champs dans un fichier MED afin d'ordonner les cellules
par type gĂŠomĂŠtrique.
L'idĂŠe est d'utiliser les deux mĂŠthodes ``MEDCouplingUMesh.sortCellsInMEDFileFrmt()`` et
``DataArrayDouble.renumberInPlace()`` pour renumĂŠroter manuellement une *copie* de ``fPart1``: ::
fPart1Cpy = fPart1.deepCopy()
o2n = fPart1Cpy.getMesh().sortCellsInMEDFileFrmt()
fPart1Cpy.getArray().renumberInPlace(o2n)
``fPart1Cpy`` est dĂŠsormais normalisĂŠ pour ĂŞtre stockĂŠ dans un fichier MED (ce que nous verrons plus loin)
VÊrifier que ``fPart1Cpy`` et ``fPart1`` sont les mêmes à une permutation près (``MEDCouplingFieldDouble.substractInPlaceDM()``) ::
fPart1Cpy.substractInPlaceDM(fPart1,12,1e-12)
fPart1Cpy.getArray().abs()
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
de ``fPart1`` et ``fPart1Cpy`` sont ĂŠgaux Ă une permutation de cellule
et/ou noeuds.
AgrĂŠger des champs
~~~~~~~~~~~~~~~~~~
AgrĂŠger ``fPart1`` et ``fPart2`` (utiliser ``MEDCouplingFieldDouble.MergeFields()``). Et mettre le rĂŠsultat de l'agrĂŠgation
dans ``fPart12``. ::
fPart12 = mc.MEDCouplingFieldDouble.MergeFields([fPart1,fPart2])
fPart12.writeVTK("ExoField_fPart12.vtu")
.. note:: La mĂŠthode ``MEDCouplingFieldDouble.MergeFields()`` devrait vraiment se
nommer ``MEDCouplingFieldDouble.AggregateFields()`` ...
.. image:: images/FieldDouble1_2.png
Evaluation d'un champ en des points donnĂŠs de l'espace
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Evaluer la valeur du champ ``fPart12`` calculĂŠ prĂŠcĂŠdemment sur les barycentres des cellules de son
maillage (variable ``bary``) et mettre le rĂŠsultat dans ``arr1``.
Utiliser pour cela les mĂŠthodes ``MEDCouplingFieldDouble.getValueOnMulti()`` et ``MEDCouplingMesh.computeCellCenterOfMass()``.
De manière similaire, Êvaluer ensuite directement le champ ``f`` en utilisant la même liste de points
que prĂŠcĂŠdemment (``bary``) et mettre le rĂŠsultat dans ``arr2``.
VĂŠrifier ensuite que ``arr1`` et ``arr2`` sont bien ĂŠgaux: ::
bary = fPart12.getMesh().computeCellCenterOfMass()
arr1 = fPart12.getValueOnMulti(bary)
arr2 = f.getValueOnMulti(bary)
delta = arr1-arr2
delta.abs()
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ĂŠ.
Pour les champs aux noeuds (P1), les cellules doivent être de types simples (triangles, tÊtraèdres) et une interpolation
linĂŠaire est alors utilisĂŠe.
.. note:: Cette technique peut ĂŞtre utilisĂŠe pour juger rapidement de la qualitĂŠ d'une interpolation.
OpĂŠrations sur les champs
~~~~~~~~~~~~~~~~~~~~~~~~~
Calculer l'intÊgrale du champ ``fPart12`` sur le maillage, et la retrouver d'une autre manière en utilisant
la mĂŠthode ``DataArrayDouble.accumulate()`` sur le tableau de valeurs de ce champ.
On rappelle que, vu le maillage simplifiĂŠ en jeu, les cellules ont toutes un volume unitĂŠ. ::
integ1 = fPart12.integral(0,True)
integ1_bis = fPart12.getArray().accumulate()[0]
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 )
Exploser un champ - Vecteurs de dĂŠplacement
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Nous allons maintenant crĂŠer un nouveau maillage reprĂŠsentant l'*ĂŠclatĂŠ* du maillage initial.
Partant du maillage ``mesh`` crĂŠer un champ vectoriel aux cellules ``fVec`` ayant 3 composantes reprĂŠsentant
le vecteur dĂŠplacement entre le point [5.,5.,5.] et le barycentre de chaque cellule du maillage.
Utiliser la mĂŠthode ``MEDCouplingMesh.fillFromAnalytic()`` : ::
fVec = mesh.fillFromAnalytic(mc.ON_CELLS,3,"(x-5.)*IVec+(y-5.)*JVec+(z-5.)*KVec")
.. note:: Les identifiants spÊciaux ``IVec``, ``JVec`` et ``KVec`` reprÊsentent les vecteurs unitaires du repère.
CrĂŠer ensuite une rĂŠduction de ``fVec`` (nommĂŠe ``fVecPart1``) sur les cellules ``ids1`` prĂŠcĂŠdemment obtenues : ::
fVecPart1 = fVec.buildSubPart(ids1)
fVecPart1.setName("fVecPart1")
Construire le champ scalaire ``fPart1Exploded`` ayant les mĂŞmes valeurs que ``fPart1`` mais reposant sur un maillage *eclatĂŠ*
par rapport Ă celui de ``fPart1.getMesh()``. Pour exploser ``fPart1.getMesh()`` utiliser le champ de dĂŠplacement vectoriel
``fVecPart1`` afin d'appliquer Ă chaque cellule la translation associĂŠe. ::
cells = fPart1.getMesh().getNumberOfCells() * [None]
for icell,vec in enumerate(fVecPart1.getArray()):
m = fPart1.getMesh()[[icell]]
m.zipCoords() # Not mandatory but saves memory
m.translate(vec)
cells[icell] = m
pass
meshFVecPart1Exploded = mc.MEDCouplingUMesh.MergeUMeshes(cells)
fPart1.setMesh(meshFVecPart1Exploded)
fPart1.writeVTK("ExoField_fPart1_explo.vtu")
Et voilĂ ce que vous devriez obtenir:
.. image:: images/FieldDouble1_1_exploded.png
:scale: 120
Solution
~~~~~~~~
:ref:`python_testMEDCouplingfielddouble1_solution`