doc/tutorial/medcoupling_3Dcube.rst

   1
   2 Meshing a 3D cube
   3 -----------------
   4
   5 Objective
   6 ~~~~~~~~~
   7
   8 The meshing class of the SALOME MED module allows user to create a mesh from scratch.
   9 In this example we propose to build a mesh on a 3D cube by three methods (classical method, method by extrusion and grid method). Each cell of the mesh must be a hexaedron.
  10 We see also how creating a group.
  11 Then we create a field on all the 3D cube.
  12 Each result will be saved in a med file, and then visualized with the SMESH module of Salome.
  13
  14 In spite of a mesh in MEDCoupling has only one dimension, it's possible to save some meshes with different dimension in one med file. We will see this method.
  15
  16 .. image:: images/Mesh_cube3D.jpg
  17
  18
  19 Beginning of implementation
  20 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
  21
  22 To implement this exercice we use the python language script and import the MEDCoupling and MEDLoader parts of the MED module. We need also mathematical functions, so we import the python math module::
  23
  24         from MEDCoupling import *
  25         from MEDLoader import *
  26         from math import *
  27
  28 You must define 3 variables for space dimension, number of nodes on each dimension and total number of nodes::
  29
  30         MeshDim3D = 3
  31         MeshDim2D = 2
  32         N = ...
  33         NbCell2D = (N-1)*(N-1)
  34         NbCell3D = NbCell2D*(N-1)
  35         NbNode2D = N*N
  36         NbNode3D = NbNode2D*N
  37
  38
  39 Classical method
  40 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  41
  42 First instanciate a meshing object. Therefore, we need to define :
  43
  44  * its name
  45  * its dimension
  46  * its number of cells
  47
  48 .. note:: All this initialisation are necessary. If one lacks, you will have a segmentation fault!.
  49
  50 ::
  51
  52         mesh=MEDCouplingUMesh.New()
  53         mesh.setMeshDimension(3)
  54         mesh.allocateCells(...)
  55         mesh.setName("3Dcube")
  56
  57 Definition of nodes coordinates
  58 ```````````````````````````````
  59
  60 Define the coordinates of the nodes of the 3D cube mesh, and then use the setCoords function to set it::
  61
  62         # Initialisation of coordinates
  63         coordinates = []
  64         for k in range(N):
  65                 for j in range(N):
  66                         for i in range(N):
  67                                 coordinates.append(...)
  68
  69         myCoords = DataArrayDouble.New()
  70         myCoords.setValues(coordinates,nbOfNodes,3)
  71         mesh.setCoords(myCoords)
  72
  73
  74 Definition of hexahedrons connectivity
  75 ``````````````````````````````````````
  76 For each hexahedron of the mesh, you have to give its connectivity: the list of the nodes which belong to the hexahedron. The order of the nodes in the connectivity array must respect the MEDCoupling format (see the following figure).
  77
  78 .. image:: images/cube.jpg
  79
  80 .. warning:: Connectivity elements begin to 0 from (n-1) elements
  81
  82 ::
  83
  84         connectivity = []
  85         for k in range(N-1):
  86                 for j in range(N-1):
  87                         for i in range(N-1):
  88                                 inode = ...
  89                                 connectivity.append(inode)
  90                                 connectivity.append(...)
  91
  92         # Adding cells in meshing
  93         for i in range(nbOfCells):
  94                 mesh.insertNextCell(NORM_HEXA8,8,connectivity[8*i:8*(i+1)])
  95                 pass
  96
  97         # Finishing insertion
  98         mesh.finishInsertingCells()
  99         mesh.checkCoherency()
 100
 101 Method by extrusion
 102 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 103
 104 In order to create a extruded mesh, we need one 2D mesh and one 1D mesh which define the vector of extrusion and the number of steps.
 105
 106 Definition of 2D mesh
 107 ``````````````````````````````````````
 108 ::
 109
 110         coordinates = []
 111         for j in range(N):
 112                 for i in range(N):
 113                         coordinates.append(...)
 114                         ...
 115         Connectivities = [...]
 116         myCoords = DataArrayDouble.New()
 117         myCoords.setValues(coordinates,NbNode2D,MeshDim2D)
 118
 119         m1 = MEDCouplingUMesh.New()
 120         m1.setMeshDimension(MeshDim2D)
 121         m1.allocateCells(NbCell2D)
 122         m1.setCoords(myCoords)
 123         m1.setName("2D_Support")
 124
 125         for i in range(NbCell2D):
 126                 m1.insertNextCell(NORM_QUAD4,4,Connectivities[4*i:4*(i+1)])
 127         m1.finishInsertingCells()
 128         m1.changeSpaceDimension(3)
 129
 130 Definition of 1D mesh
 131 ``````````````````````````````````````
 132 ::
 133
 134         coords = [ ... ]
 135         conn   = [ ... ]
 136         m2 = MEDCouplingUMesh.New()
 137         m2.setMeshDimension(1)
 138         m2.allocateCells(3)
 139         m2.insertNextCell(NORM_SEG2,2,conn[0:2])
 140         m2.insertNextCell(NORM_SEG2,2,conn[2:4])
 141         m2.insertNextCell(NORM_SEG2,2,conn[4:6])
 142         m2.finishInsertingCells()
 143         myCoords1D=DataArrayDouble.New()
 144         myCoords1D.setValues(coords,4,1)
 145         m2.setCoords(myCoords1D)
 146         m2.changeSpaceDimension(3)
 147
 148 Definition of extruded mesh
 149 ``````````````````````````````````````
 150
 151 Since 1D meshing isn't well oriented (along 0x vector), we need to imply a transformation on it.
 152 Then, we can extrude 2D meshing.
 153
 154 ::
 155
 156         center = [...]
 157         vector = [...]
 158         m2.rotate(...)
 159         m3 = m1.buildExtrudedMesh(m2,0)
 160         m3.setName("Extrusion")
 161
 162 Grid method
 163 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 164
 165 it's the easiest way to create a grid since you have no connectivity to set. They will be automatically setting. Incrementation of ids will be made first along Ox axis, then along Oy axis and finally along Oz axis.
 166
 167 ::
 168
 169         mesh=MEDCouplingCMesh.New()
 170         coordsX=DataArrayDouble.New()
 171         arrX=[ ... ]
 172         coordsX.setValues(arrX,4,1)
 173         coordsY=DataArrayDouble.New()
 174         arrY=[ ... ]
 175         coordsY.setValues(arrY,4,1)
 176         coordsZ=DataArrayDouble.New()
 177         arrZ=[ ... ]
 178         coordsZ.setValues(arrZ,4,1)
 179         mesh.setCoords(coordsX,coordsY,coordsZ)
 180
 181 Really in order to save this mesh, you need to transform this structured mesh to an unstructerd mesh.
 182 ::
 183
 184         meshU=mesh.buildUnstructured()
 185         meshU.setName("Grid")
 186
 187 Creation of groups
 188 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 189
 190 A group is a set of cells defining by their id. This set must the input for creating a group.
 191 Generally ids cells using in group are known. So you just need put these ids in a DataArray.
 192 ::
 193
 194         tabIdCells = DataArrayInt.New()
 195         IdCells = [ ... ]
 196         tabIdCells.setValues(IdCells,...)
 197
 198
 199 .. note:: It's also possible to retrieve ids cells from a submesh of the principal mesh.
 200
 201 ::
 202
 203         ret,tabIdCells = mesh.areCellsIncludedIn(subMesh,0)
 204
 205
 206 Once the DataArray is created, some initializations have to be done.
 207 ::
 208
 209         # Definition of the name group
 210         tabIdCells.setName("meshGroup")
 211
 212
 213 In order to add a group on a mesh, you need to transform your unstructured mesh in a file unstructured mesh.
 214 Moreover, we need to define:
 215
 216  * its name
 217  * its description
 218  * its coordinates
 219  * its dimension
 220  * its number of cells
 221
 222
 223 ::
 224
 225         # Passing MEDCoupling to MEDFile
 226         fmeshU = MEDFileUMesh.New()
 227         fmeshU.setName("Grid")
 228         fmeshU.setDescription("IHopeToConvinceLastMEDMEMUsers")
 229         myCoords = meshU.getCoords()
 230         fmeshU.setCoords(myCoords)
 231         fmeshU.setMeshAtLevel(0,meshU)
 232
 233 Then, you can
 234 Therefore, you need to define the level (ie. the dimension) of the group.
 235 This information is given by a number : 0,-1 or -2.
 236
 237  * 0 means the same level at mesh
 238
 239 ::
 240
 241         fmeshU.setGroupsAtLevel(0,[tabIdCells],False)
 242
 243 Create field on 3D cube
 244 ~~~~~~~~~~~~~~~~~~~~~~~
 245
 246 First you need to create a CouplingField and initialize some data:
 247
 248  * its name
 249  * its support (ie mesh)
 250  * its nature
 251  * its values
 252
 253
 254 The field will be a sin function dependant of distance of the barycenter of each cell from origin. So we need to create a barycenter field on the 3D mesh::
 255
 256         # Creation of field : with following definition
 257         # => Definition of the mesh support
 258         # => Definition of field name
 259         # => Definition of field nature
 260         field = MEDCouplingFieldDouble.New(ON_CELLS)
 261         field.setMesh(mesh)
 262         field.setName("field")
 263         field.setNature(Integral)
 264
 265         # Computing and setting field values
 266         myCoords=DataArrayDouble.New()
 267         sampleTab=[]
 268         bar = mesh.getBarycenterAndOwner()
 269         print bar.getNbOfElems()
 270         for i in range(nbOfCells):
 271                 x = bar.getIJ(...)
 272                 y = bar.getIJ(...)
 273                 z = bar.getIJ(...)
 274                 d = sqrt(x*x+y*y+z*z)
 275                 sinus = sin(d)
 276         .       sampleTab.append(sinus)
 277
 278         myCoords.setValues(sampleTab,nbOfCells,1)
 279         field.setArray(myCoords)
 280
 281
 282 Saving the mesh in a med file
 283 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 284 General Case
 285 ````````````
 286 ::
 287
 288         medFileName = "MEDCoupling_Extrudedcube3D.med"
 289         MEDLoader.WriteUMesh(medFileName,meshU,True)
 290
 291 .. note:: True / False in Write* functions : True for overwriting existing file and False for adding in existing file
 292
 293 Multi mesh Case
 294 ````````````````
 295
 296 In spite of a MEDCoupling mesh has only one dimension, it's possible to genrate a file with multi dimension.
 297 Therefore, you need to create as meshes as necessary dimensions.
 298
 299 You have to give the connectivity of the faces on the bottom face of the 3D cube: the list of the nodes which belong to the face.
 300 The connectivity must respect following figure:
 301
 302 .. image:: images/face.jpg
 303
 304 ::
 305
 306         # Extraction of surfacic meshing
 307         pt=[0.,0.,0.]
 308         vec=[0.,0.,1.]
 309         nodes = mesh.findNodesOnPlane(pt,vec,1e-12)
 310         mesh2D = mesh.buildFacePartOfMySelfNode(nodes,True)
 311         #print mesh2D
 312         mesh2D.setName("3Dcube")
 313         mesh2D.checkCoherency()
 314
 315         medFileName = "MEDCoupling_cube3D.med"
 316         meshes=[mesh2D,mesh]
 317         MEDLoader.WriteUMeshes(medFileName,meshes,True)
 318
 319 Group Case
 320 ````````````
 321 ::
 322
 323         medFileName = "MEDCoupling_Gridcube3D.med"
 324         fmeshU.write(medFileName,2)
 325
 326 Saving the fields in the med file
 327 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 328
 329 ::
 330
 331         MEDLoader.WriteField(medFileName,field,False)
 332
 333 Visualize the mesh with the SMESH module of Salome
 334 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 335
 336 Launch Salome platform, then select SMESH module and import your MED file. First You can see the number of elements of your mesh. For that, select your mesh in the object browser, set select Mesh menu and "Advanced Mesh Info" option. Verify the number of faces and the number of hexahedrons. To visualize your mesh: click right bottom on your mesh and select "Show" option. You can also visualize your groups. Select one group, click right bottom on your group and select "Show only" option.
 337
 338 Visualize the fields with the VISU module of Salome
 339 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 340
 341 Launch Salome platform, then select VISU module and import your MED file. You can see in the object browser the 2 fields you have created. Then you have to create a scalar map on each field to visualize them.
 342
 343 .. image:: images/Field_Cube3D.jpg
 344
 345 Solution
 346 ~~~~~~~~
 347
 348 :ref:`python_testMEDCouplingcube_solution`