src/ElectromagnetismFluxDisc/plugin/Test/test_flux_disc0.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 #  -*- coding: utf-8 -*-
  21
  22 from paraview.simple import *
  23 from vtk.util import numpy_support
  24 import numpy as np
  25 import vtk
  26
  27 from medcoupling import *
  28
  29 def MyAssert(clue):
  30     if not clue:
  31         raise RuntimeError("Assertion failed !")
  32
  33 def Write(ds,fn):
  34     writer = vtk.vtkUnstructuredGridWriter()
  35     writer.SetInputData(ds)
  36     writer.SetFileName(fn)
  37     writer.Update()
  38
  39 def MCField(fn,fieldName):
  40     mm = MEDFileMesh.New(fn)
  41     fs = MEDFileFields.New(fn)
  42     return fs[fieldName][0].field(mm)
  43
  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]
  51
  52 ref_flux_b_a = -60.47
  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]
  57
  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']
  61 mr.UpdatePipeline()
  62
  63 resRadial = 200
  64 resCircum = 200
  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
  72
  73 mw = MEDWriter(FileName="Base.med",Input=fluxDisc1)
  74 mw.UpdatePipeline()
  75
  76 ds0 = servermanager.Fetch(fluxDisc1)
  77 MyAssert(ds0.GetNumberOfCells()==resRadial*resCircum)
  78 MyAssert(ds0.GetNumberOfPoints()==resRadial*resCircum+1)
  79
  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))
  82
  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))
  85
  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))
  88
  89 Write(ds0.GetBlock(0),"Base.vtu")
  90 mc_f_base = MCField("Base.med","H_HsOmega")
  91
  92 # On inverse l'axe et on verifie que le flux est inversé
  93
  94 fluxDisc1.SliceType.Axis = [-elt for elt in vector]
  95 mw = MEDWriter(FileName="Invert.med",Input=fluxDisc1)
  96 mw.UpdatePipeline()
  97
  98 ds1 = servermanager.Fetch(fluxDisc1)
  99
 100 MyAssert(ds1.GetNumberOfCells()==resRadial*resCircum)
 101 MyAssert(ds1.GetNumberOfPoints()==resRadial*resCircum+1)
 102
 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))
 105
 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))
 108
 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))
 111
 112
 113 Write(ds0.GetBlock(0),"Invert.vtu")
 114 mc_f_invert = MCField("Invert.med","H_HsOmega")
 115
 116 ###
 117
 118 ComputeFlux(mc_f_base)
 119 ComputeFlux(mc_f_invert)
 120
 121 # debug : pourquoi une telle difference de flux sur H_HsOmega ?
 122
 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))
 126
 127 # OK pour l'area
 128
 129 m_base = mc_f_base.getMesh()
 130 m_invert = mc_f_invert.getMesh()
 131 base_base = DataArrayDouble.GiveBaseForPlane(vector)
 132
 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
 136
 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])
 139
 140 m_base.setCoords(base_new_coords)
 141 m_invert.setCoords(invert_new_coords)
 142
 143 mc_f_base.write("base_dbg.vtu")
 144 mc_f_invert.write("invert_dbg.vtu")
 145
 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)))
 150
 151 a2 = a[:] ; a2.sort() ; assert( a2.isIota(len(a2) ) )
 152
 153 array_base_in_invert_ref = mc_f_base.getArray()[a]
 154
 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.)
 157
 158 ze_ids = ok_ids
 159 ze_ids = not_of_ids
 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]))
 163
 164 mc_f_invert[ze_ids].write("invert_dbg_pb.vtu")