src/RateOfFlowThroughSection/script/calcul_sediment_deposit.py

   1 # Copyright (C) 2021  CEA/DEN, EDF R&D
   2 #
   3 # This library is free software; you can redistribute it and/or
   4 # modify it under the terms of the GNU Lesser General Public
   5 # License as published by the Free Software Foundation; either
   6 # version 2.1 of the License, or (at your option) any later version.
   7 #
   8 # This library is distributed in the hope that it will be useful,
   9 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 # Lesser General Public License for more details.
  12 #
  13 # You should have received a copy of the GNU Lesser General Public
  14 # License along with this library; if not, write to the Free Software
  15 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 #
  17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 #
  19
  20 from medcoupling import *
  21
  22 data_file="Res_sisy.med"
  23 poly_file="contour.med"
  24 cutoff=0.
  25
  26 polygon=ReadMeshFromFile(poly_file)
  27 data=MEDFileMesh.New(data_file)
  28 mesh=data[0]
  29 mesh.changeSpaceDimension(2,0.)
  30 mesh=mesh.deepCopy()
  31 assert(polygon.getCoords()[:,2].isUniform(0,1e-12))
  32 polygon.changeSpaceDimension(2,0.)
  33 polygon.mergeNodes(1e-12) #
  34 polygon_1sgt=MEDCoupling1SGTUMesh(polygon)
  35 conn=polygon_1sgt.getNodalConnectivity().deepCopy()
  36 conn.rearrange(2)
  37 notNullCellsPolygon=(conn[:,0]-conn[:,1]).findIdsNotEqual(0)
  38 polygon=polygon_1sgt.buildUnstructured()[notNullCellsPolygon]
  39 polygon_2d=MEDCouplingUMesh("mesh",2)
  40 polygon_2d.setCoords(polygon.getCoords().deepCopy())
  41 polygon_2d.allocateCells()
  42 conn=MEDCoupling1SGTUMesh(polygon).getNodalConnectivity()
  43 conn.rearrange(2)
  44 conn2=conn.fromLinkedListOfPairToList()
  45 assert(conn2[0]==conn2[-1])
  46 conn2.popBackSilent()
  47 polygon_2d.insertNextCell(NORM_POLYGON,conn2.getValues())
  48 clockWise = polygon_2d.getMeasureField(False).getIJ(0,0) > 0.
  49 #
  50 side={True : 1 , False : 0}[clockWise]
  51
  52 (Xmin,Xmax),(Ymin,Ymax)=mesh.getBoundingBox()
  53 TmpCenter=( (Xmin+Xmax)/2., (Ymin+Ymax)/2. )
  54 Tmpalpha=1/max(Xmax-Xmin,Ymax-Ymin)
  55 mesh.translate(-DataArrayDouble(TmpCenter,1,2))
  56 mesh.scale([0.,0.],Tmpalpha)
  57 polygon.translate(-DataArrayDouble(TmpCenter,1,2))
  58 polygon.scale([0.,0.],Tmpalpha)
  59 mesh2,line_inter,cellid_in_2d,cellid_in1d = MEDCouplingUMesh.Intersect2DMeshWith1DLine(mesh,polygon,1e-12)
  60 coo=mesh2.getCoords().deepCopy()
  61 coo2=coo*(1/Tmpalpha)+TmpCenter
  62 mesh2.setCoords(coo2)
  63
  64 twodcells_to_remove = cellid_in1d[:,side]
  65 twodcells_to_keep = cellid_in1d[:,(side+1)%2]
  66 ids=twodcells_to_keep.findIdsNotEqual(-1)
  67 twodcells_to_keep=twodcells_to_keep[ids]
  68 hotspot = twodcells_to_keep[0]
  69
  70 twodcells_to_keep = twodcells_to_keep[ids] # les cells2D de bord du domaine dans le referentiel de sortie 2D
  71 twodcells_to_remove.sort()
  72 ids=twodcells_to_remove.findIdsNotEqual(-1)
  73 twodcells_to_remove=twodcells_to_remove[ids]
  74 allcells=twodcells_to_remove.buildComplement(mesh2.getNumberOfCells())
  75 mesh2_without_cells_around_polygon=mesh2[allcells]
  76 grps=mesh2_without_cells_around_polygon.partitionBySpreadZone()
  77 grps=[allcells[elt] for elt in grps]
  78 assert(len(grps)==2)
  79 zeGrp = None
  80 if (hotspot in grps[0]) and (hotspot not in grps[1]):
  81     zeGrp = grps[0]
  82 if (hotspot not in grps[0]) and (hotspot in grps[1]):
  83     zeGrp = grps[1]
  84 if not zeGrp:
  85     raise RuntimeError("Ooops")
  86     pass
  87 mesh3 = mesh2[zeGrp]
  88 totVol = mesh3.getMeasureField(True).accumulate()[0]
  89 refVol = polygon_2d.getMeasureField(True).accumulate()[0]
  90 assert(abs((totVol-refVol)/refVol)<1e-6)
  91 #
  92 original_cell_ids_2d=cellid_in_2d[zeGrp] ; original_cell_ids_2d.sort() # les cells dans le referentiel original 2D
  93 all_cut_2d_cells=cellid_in1d[:]
  94 all_cut_2d_cells.rearrange(1)
  95 all_cut_2d_cells.sort()
  96 all_cut_2d_cells=all_cut_2d_cells.buildUnique() # les cells qui ont subit un split dans le referentiel de sortie 2D
  97
  98 all_cut_2d_cells_origin=cellid_in_2d[all_cut_2d_cells]
  99 untouched_2d_cells=original_cell_ids_2d.buildSubstraction(all_cut_2d_cells_origin)
 100
 101 cells_at_boundary=all_cut_2d_cells.buildIntersection(zeGrp) # les cellules decoupées dans le referentiel de sortie 2D
 102 cells_at_boundary_origin=cellid_in_2d[cells_at_boundary]
 103 mesh3=mesh2[cells_at_boundary]
 104 vol = mesh3.getMeasureField(True).getArray()
 105 #volRef = mesh[cells_at_boundary_origin].getMeasureField(True).getArray()
 106 #centers=mesh3.computeCellCenterOfMass()
 107 #
 108 tmp0=mesh[untouched_2d_cells]
 109 (Xmin,Xmax),(Ymin,Ymax)=mesh3.getBoundingBox()
 110 ZeCenter=( (Xmin+Xmax)/2., (Ymin+Ymax)/2. )
 111 alpha=1/max(Xmax-Xmin,Ymax-Ymin)
 112 mesh3.translate(-DataArrayDouble(ZeCenter,1,2))
 113 mesh3.scale([0.,0.],alpha)
 114 centers_around_zero=mesh3.computeCellCenterOfMass()
 115 centers=centers_around_zero*(1/alpha)+ZeCenter
 116 #
 117 evolution_multiTS=MEDFileFloatFieldMultiTS(data_file,"EVOLUTION",False)
 118 for ts in range(10):
 119     evolution_1TS=evolution_multiTS[ts]
 120     evolution_1TS.loadArrays()
 121     f=evolution_1TS.field(data).convertToDblField()
 122     f_easy=f[untouched_2d_cells]
 123     f_easy.write("f_easy.med") #
 124     weights = f_easy.getDiscretization().getMeasureField(f_easy.getMesh(),True).getArray()
 125     positive_ids = f_easy.getArray().findIdsGreaterThan(cutoff)
 126     negative_ids = positive_ids.buildComplement(f_easy.getMesh().getNumberOfNodes())
 127     positive_part = (f_easy.getArray()[positive_ids]*weights[positive_ids]).accumulate(0)
 128     negative_part = (f_easy.getArray()[negative_ids]*weights[negative_ids]).accumulate(0)
 129     #
 130     f_hard = f[cells_at_boundary_origin]
 131     hard_part = f_hard.getValueOnMulti(centers)
 132     positive_ids_hard = hard_part.findIdsGreaterThan(cutoff)
 133     negative_ids_hard = positive_ids_hard.buildComplement(len(cells_at_boundary_origin))
 134     positive_part_hard=(hard_part[positive_ids_hard]*vol[positive_ids_hard]).accumulate(0)
 135     negative_part_hard=(hard_part[negative_ids_hard]*vol[negative_ids_hard]).accumulate(0)
 136     print("Time step %d (%ld)"%(ts,evolution_1TS.getTime()[-1]),positive_part+positive_part_hard+negative_part+negative_part_hard)
 137     evolution_1TS.unloadArrays()
 138     pass