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 # -*- coding: utf-8 -*-
22 from paraview.simple import *
23 from vtk.util import numpy_support
27 from medcoupling import *
31 raise RuntimeError("Assertion failed !")
34 writer = vtk.vtkUnstructuredGridWriter()
35 writer.SetInputData(ds)
36 writer.SetFileName(fn)
39 def MCField(fn,fieldName):
40 mm = MEDFileMesh.New(fn)
41 fs = MEDFileFields.New(fn)
42 return fs[fieldName][0].field(mm)
44 def ComputeFlux(mc_f):
45 area = mc_f.getMesh().getMeasureField(True).getArray()
46 mc_vector = DataArrayDouble(vector,1,3)
47 mc_vector /= mc_vector.magnitude()[0]
48 mc_vector = DataArrayDouble.Aggregate(mc_f.getNumberOfTuples()*[mc_vector])
49 flux_per_cell = DataArrayDouble.Dot(mc_f.getArray(),mc_vector)*area
50 return flux_per_cell.accumulate()[0]
53 ref_flux_h_a = -127446
54 ref_flux_hsomega = 1896.991081794051
55 vector = [0.13,-0.47,0.87]
56 center = [0.0, 0.0, 3900.0]
58 mr = MEDReader(FileName='/home/H87074/TMP81/paravistests_new/FilesForTests/Electromagnetism/mesh_benjamin_8_sept_2020.med')
59 mr.AllArrays = ['TS0/Mesh_1/ComSup0/B_A@@][@@P0', 'TS0/Mesh_1/ComSup0/B_HsOmega@@][@@P0', 'TS0/Mesh_1/ComSup0/H_A@@][@@P0', 'TS0/Mesh_1/ComSup0/H_HsOmega@@][@@P0']
60 mr.AllTimeSteps = ['0000', '0001', '0002', '0003', '0004', '0005', '0006', '0007', '0008', '0009']
65 fluxDisc1 = FluxDisc(Input=mr)
66 fluxDisc1.SliceType = 'Cylinder'
67 fluxDisc1.SliceType.Center = center
68 fluxDisc1.SliceType.Axis = vector
69 fluxDisc1.SliceType.Radius = 293
70 fluxDisc1.RadialResolution = resRadial
71 fluxDisc1.CircumferentialResolution = resCircum
73 mw = MEDWriter(FileName="Base.med",Input=fluxDisc1)
76 ds0 = servermanager.Fetch(fluxDisc1)
77 MyAssert(ds0.GetNumberOfCells()==resRadial*resCircum)
78 MyAssert(ds0.GetNumberOfPoints()==resRadial*resCircum+1)
80 flux_b_a = ds0.GetBlock(0).GetCellData().GetArray("B_A_flux").GetValue(0)
81 MyAssert(np.isclose(ref_flux_b_a,flux_b_a,0.1,0))
83 flux_h_a = ds0.GetBlock(0).GetCellData().GetArray("H_A_flux").GetValue(0)
84 #MyAssert(np.isclose(ref_flux_h_a,flux_h_a,0.01,0))
86 flux_hsomega = ds0.GetBlock(0).GetCellData().GetArray("H_HsOmega_flux").GetValue(0)
87 #MyAssert(np.isclose(ref_flux_hsomega,flux_hsomega,0.01,0))
89 Write(ds0.GetBlock(0),"Base.vtu")
90 mc_f_base = MCField("Base.med","H_HsOmega")
92 # On inverse l'axe et on verifie que le flux est inversé
94 fluxDisc1.SliceType.Axis = [-elt for elt in vector]
95 mw = MEDWriter(FileName="Invert.med",Input=fluxDisc1)
98 ds1 = servermanager.Fetch(fluxDisc1)
100 MyAssert(ds1.GetNumberOfCells()==resRadial*resCircum)
101 MyAssert(ds1.GetNumberOfPoints()==resRadial*resCircum+1)
103 flux_b_a = ds1.GetBlock(0).GetCellData().GetArray("B_A_flux").GetValue(0)
104 #MyAssert(np.isclose(-ref_flux_b_a,flux_b_a,0.1,0))
106 flux_h_a = ds1.GetBlock(0).GetCellData().GetArray("H_A_flux").GetValue(0)
107 #MyAssert(np.isclose(-ref_flux_h_a,flux_h_a,0.1,0))
109 flux_hsomega = ds1.GetBlock(0).GetCellData().GetArray("H_HsOmega_flux").GetValue(0)
110 #MyAssert(np.isclose(-ref_flux_hsomega,flux_hsomega,0.1,0))
113 Write(ds0.GetBlock(0),"Invert.vtu")
114 mc_f_invert = MCField("Invert.med","H_HsOmega")
118 ComputeFlux(mc_f_base)
119 ComputeFlux(mc_f_invert)
121 # debug : pourquoi une telle difference de flux sur H_HsOmega ?
123 a = mc_f_base.getMesh().getMeasureField(True).getArray()
124 b = mc_f_invert.getMesh().getMeasureField(True).getArray()
125 assert(a.isEqual(b,1e-10))
129 m_base = mc_f_base.getMesh()
130 m_invert = mc_f_invert.getMesh()
131 base_base = DataArrayDouble.GiveBaseForPlane(vector)
133 newX = DataArrayDouble(resRadial*resCircum+1,3) ; newX[:] = base_base[0]
134 newY = DataArrayDouble(resRadial*resCircum+1,3) ; newY[:] = base_base[1]
135 newZ = DataArrayDouble(resRadial*resCircum+1,3) ; newZ[:] = 0
137 base_new_coords=DataArrayDouble.Meld([DataArrayDouble.Dot(m_base.getCoords(),newX),DataArrayDouble.Dot(m_base.getCoords(),newY),newZ])
138 invert_new_coords=DataArrayDouble.Meld([DataArrayDouble.Dot(m_invert.getCoords(),newX),DataArrayDouble.Dot(m_invert.getCoords(),newY),newZ])
140 m_base.setCoords(base_new_coords)
141 m_invert.setCoords(invert_new_coords)
143 mc_f_base.write("base_dbg.vtu")
144 mc_f_invert.write("invert_dbg.vtu")
146 m_base.changeSpaceDimension(2,0.)
147 m_invert.changeSpaceDimension(2,0.)
148 a,b= m_base.getCellsContainingPoints(m_invert.computeCellCenterOfMass(),1e-12)
149 assert(b.isIota(len(b)))
151 a2 = a[:] ; a2.sort() ; assert( a2.isIota(len(a2) ) )
153 array_base_in_invert_ref = mc_f_base.getArray()[a]
155 not_of_ids = (array_base_in_invert_ref-mc_f_invert.getArray()).magnitude().findIdsGreaterThan(1.)
156 ok_ids = (array_base_in_invert_ref-mc_f_invert.getArray()).magnitude().findIdsLowerThan(1.)
160 print(ComputeFlux(mc_f_invert[ze_ids]))
161 ze_ids_base_ref = a.findIdForEach(ze_ids)
162 print(ComputeFlux(mc_f_base[ze_ids_base_ref]))
164 mc_f_invert[ze_ids].write("invert_dbg_pb.vtu")