CDMATH/tests/ressources/2DBrickWall/create_mesh_brickWall.py

   1 # -*-coding:utf-8 -*
   2
   3 #### Maillage d'un rectangle par une grille décalée (mur de briques)
   4 ### input : xmin, xmax, nx, ymin, ymax, ny
   5 ### output : squareWithBrickWall.vtu squareWithBrickWall.med
   6
   7 import medcoupling as mc
   8 import math
   9 import MEDLoader as ML
  10
  11 def createBrickWallMesh( xmin=0., xmax=1., nx=15, ymin=0., ymax=1., ny=15,mesh_name="squareWithBrickWall"):
  12
  13     dx = (xmax-xmin)/nx
  14     dy=(ymax-ymin)/ny
  15
  16     print( "Creating BrickWall mesh with nx=",nx,"ny=",ny, "nb cells=",nx*ny )
  17
  18     # Building the initial rectangular cell, centered at 0,0
  19     d = mc.DataArrayDouble(4,2)
  20     d[0,0] = -dx/2
  21     d[0,1] =  dy/2
  22     d[1,0] =  dx/2
  23     d[1,1] =  dy/2
  24     d[2,0] =  dx/2
  25     d[2,1] = -dy/2
  26     d[3,0] = -dx/2
  27     d[3,1] = -dy/2
  28
  29     d.setInfoOnComponents(["X [m]","Y [m]"])
  30
  31     print( "Uniform array ?", d.magnitude().isUniform(0.5*math.sqrt(dx*dx+dy*dy),1e-10) )
  32
  33     # translation of the first cell
  34     translationToPerform = [[(0.5*(1+j%2)+i)*dx,(0.5+j)*dy] for i in range(nx) for j in range(ny)]
  35
  36     ds = len(translationToPerform)*[None]
  37     for pos,t in enumerate(translationToPerform):
  38                      ds[pos] = d[:]         # Perform a deep copy of d and place it at position 'pos' in ds
  39                      ds[pos] += t             # Adding a vector to a set of coordinates does a translation
  40                      pass
  41
  42     d2 = mc.DataArrayDouble.Aggregate(ds)
  43     # Build an unstructured mesh representing the final pattern
  44     mesh = mc.MEDCouplingUMesh(mesh_name,2)
  45     mesh.setCoords(d2)
  46     print( "Mesh dimension is", mesh.getMeshDimension() )
  47     print( "Spatial dimension is", mesh.getCoords().getNumberOfComponents() )
  48     mesh.allocateCells(nx*ny)
  49     for i in range(nx*ny):
  50             cell_connec = [4*i,4*i+1,4*i+2,4*i+3]
  51             mesh.insertNextCell(mc.NORM_QUAD4, cell_connec)
  52             pass
  53
  54     # Identifying duplicate nodes
  55     oldNbOfNodes=mesh.getNumberOfNodes()
  56     arr, areNodesMerged, newNbOfNodes=mesh.mergeNodes(1e-10)
  57     print( "oldNbOfNodes=",oldNbOfNodes,"newNbOfNodes",newNbOfNodes )
  58
  59     # Crée des polygones pour rendre conforme les mailles
  60     mesh.conformize2D(1e-10)
  61     # Check that everything is coherent (will throw if not)
  62     mesh.checkConsistencyLight()
  63
  64     # Crée les éléments 1D pour pouvoir imposer les conditions aux limites
  65     mesh_1d = mesh.computeSkin()
  66
  67     # Identifie les segments de chaque côté pour créer les groupes
  68     tol = 1e-10
  69
  70     # PB: getCellsInBoundingBox renvoie aussi les segments qui touchent la bounding box
  71     # => On boucle sur les coordonnées des barycentres
  72
  73     barycenters = mesh_1d.computeIsoBarycenterOfNodesPerCell()
  74     ids_left = []
  75     ids_right = []
  76     ids_bottom = []
  77     ids_top = []
  78     for i, coord in enumerate(barycenters):
  79         x, y = coord
  80         if abs(y-ymin) < tol :
  81           ids_bottom.append(i)
  82         elif abs(y-ymax) < tol :
  83           ids_top.append(i)
  84         elif abs(x-xmax) < tol or abs(x-xmax-dx/4) < tol or abs(x-xmax-dx/2) < tol:
  85           ids_right.append(i)
  86         elif abs(x-xmin) < tol or abs(x-xmin-dx/4) < tol or abs(x-xmin-dx/2) < tol:
  87           ids_left.append(i)
  88         else:
  89             raise ValueError("Pb with boundary construction : barycenter does not belong to any border group")
  90
  91     arr_left = mc.DataArrayInt64(ids_left)
  92     arr_right = mc.DataArrayInt64(ids_right)
  93     arr_bottom = mc.DataArrayInt64(ids_bottom)
  94     arr_top = mc.DataArrayInt64(ids_top)
  95
  96     arr_left.setName("Left")
  97     arr_right.setName("Right")
  98     arr_bottom.setName("Bottom")
  99     arr_top.setName("Top")
 100
 101     # Trie les cellules par type conformément à la convention MED fichier
 102     o2n = mesh.sortCellsInMEDFileFrmt()
 103     meshMEDFile=ML.MEDFileUMesh.New()
 104     # Ecrit le maillage 2D
 105     meshMEDFile.setMeshAtLevel(0,mesh)
 106     # Ecrit le maillage 1D
 107     meshMEDFile.setMeshAtLevel(-1,mesh_1d)
 108     # Ecrit les groupes
 109     meshMEDFile.addGroup(-1, arr_left)
 110     meshMEDFile.addGroup(-1, arr_right)
 111     meshMEDFile.addGroup(-1, arr_bottom)
 112     meshMEDFile.addGroup(-1, arr_top)
 113
 114     filename = mesh_name+str(nx*ny)+".med"
 115     # Write the result into a VTU file that can be read with ParaView
 116     #mesh.writeVTK(mesh_name+str(nx*ny)+".vtu")
 117     # Write the result into a MED file that can be read with Salomé
 118     meshMEDFile.write(filename,2) # 2 stands for write from scratch
 119
 120 if __name__ == '__main__':
 121   createBrickWallMesh(0.,1.,5 ,0.,1.,5)
 122   createBrickWallMesh(0.,1.,15,0.,1.,15)
 123   createBrickWallMesh(0.,1.,31,0.,1.,31)
 124   createBrickWallMesh(0.,1.,61,0.,1.,61)
 125   createBrickWallMesh(0.,1.,101,0.,1.,101)
 126