1 # Copyright (C) 2021 CEA/DEN, EDF R&D
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.
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.
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
17 # See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 from medcoupling import *
22 data_file="Res_sisy.med"
23 poly_file="contour.med"
26 polygon=ReadMeshFromFile(poly_file)
27 data=MEDFileMesh.New(data_file)
29 mesh.changeSpaceDimension(2,0.)
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()
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()
44 conn2=conn.fromLinkedListOfPairToList()
45 assert(conn2[0]==conn2[-1])
47 polygon_2d.insertNextCell(NORM_POLYGON,conn2.getValues())
48 clockWise = polygon_2d.getMeasureField(False).getIJ(0,0) > 0.
50 side={True : 1 , False : 0}[clockWise]
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
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]
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]
80 if (hotspot in grps[0]) and (hotspot not in grps[1]):
82 if (hotspot not in grps[0]) and (hotspot in grps[1]):
85 raise RuntimeError("Ooops")
88 totVol = mesh3.getMeasureField(True).accumulate()[0]
89 refVol = polygon_2d.getMeasureField(True).accumulate()[0]
90 assert(abs((totVol-refVol)/refVol)<1e-6)
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
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)
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()
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
117 evolution_multiTS=MEDFileFloatFieldMultiTS(data_file,"EVOLUTION",False)
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)
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()