From d151468b0a5644acb663acf7ed48b3456985ace9 Mon Sep 17 00:00:00 2001
From: Anthony Geay <anthony.geay@edf.fr>
Date: Wed, 22 Sep 2021 09:04:00 +0200
Subject: [PATCH] [EDF24091] : Use of area per cell field was missing in
 computation of Force in TorseurCIH filter

---
 src/TorseurCIH/Test.py                             | 14 +++++++++++---
 .../plugin/TorseurCIHModule/vtkTorseurCIH.cxx      |  3 +++
 2 files changed, 14 insertions(+), 3 deletions(-)

diff --git a/src/TorseurCIH/Test.py b/src/TorseurCIH/Test.py
index dfff8c4..efbf7ec 100644
--- a/src/TorseurCIH/Test.py
+++ b/src/TorseurCIH/Test.py
@@ -40,11 +40,16 @@ def fff(initialVect,normalFace,OM,area_field):
     for expo in range(5):
         pos,p,v = fffff(initialVect,normalFace,OM,area_field,[pos+i*(10**-expo) for i in range(-9,10)])
     return pos,p,v
-    
-fname = "slice.med"
+
+"""
+EDF24091 : probleme de calcul dans TorseurCIH
+"""
+
+fname = "Case_22_09_21.med"
+fieldName = "LIN_____SIEF_NOEU"
 mm=MEDFileMesh.New(fname)
 m=mm[0]
-f1ts = MEDFileField1TS(fname,"RESUME__SIEF_NOEU")
+f1ts = MEDFileField1TS(fname,fieldName)
 f = f1ts.field(mm)
 m = f.getMesh()
 area_field = m.getMeasureField(True)
@@ -69,6 +74,7 @@ F_y = matrix[:,3]*eqn[:,0] + matrix[:,1]*eqn[:,1] + matrix[:,5]*eqn[:,2]
 F_z = matrix[:,4]*eqn[:,0] + matrix[:,5]*eqn[:,1] + matrix[:,2]*eqn[:,2]
 #
 F = DataArrayDouble.Meld([F_x,F_y,F_z])
+F[:] *= area_vector
 #
 ZeForce = DataArrayDouble(F.accumulate(),1,3)
 normalFace = DataArrayDouble(eqn.accumulate(),1,3)
@@ -120,3 +126,5 @@ f.setArray(new_dist_to_base_X)
 f.setName("dist")
 f.writeVTK("distEig.vtu")"""
 #mat*w[:,0]
+
+print("{:g}".format(ForceNormale[0,0]))
diff --git a/src/TorseurCIH/plugin/TorseurCIHModule/vtkTorseurCIH.cxx b/src/TorseurCIH/plugin/TorseurCIHModule/vtkTorseurCIH.cxx
index 9ac9f26..f9f43b9 100644
--- a/src/TorseurCIH/plugin/TorseurCIHModule/vtkTorseurCIH.cxx
+++ b/src/TorseurCIH/plugin/TorseurCIHModule/vtkTorseurCIH.cxx
@@ -702,7 +702,10 @@ vtkSmartPointer<vtkTable> ComputeTorseurCIH(vtkUnstructuredGrid* usgIn)
     F_y = ForceBuilder({ 3, 1, 5 }, matrix, eqn);
     F_z = ForceBuilder({ 4, 5, 2 }, matrix, eqn);
   }
+  //
+  //
   MCAuto<DataArrayDouble> F(DataArrayDouble::Meld({ F_x, F_y, F_z }));
+  F->multiplyEqual(area_vector);
   double ZeForce[3], normalFace[3];
   F->accumulate(ZeForce);
   eqn->accumulate(normalFace);
-- 
2.30.2