From e079fe6f11ed8872b3d58ccb298252550dc0ec38 Mon Sep 17 00:00:00 2001
From: geay <anthony.geay@cea.fr>
Date: Thu, 6 Mar 2014 10:10:33 +0100
Subject: [PATCH] correct hexa27 gauss points shape functions

---
 .../GaussPoints/InterpKernelGaussCoords.cxx   | 24 +++---
 src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 27 +++++++
 src/MEDLoader/Swig/MEDLoaderTest4.py          | 79 +++++++++++++++++++
 3 files changed, 118 insertions(+), 12 deletions(-)

diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
index 97ccac266..13f90bc69 100644
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
@@ -2620,18 +2620,18 @@ void GaussInfo::hexa27aInit()
    funValue[5] =0.125*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
    funValue[6] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
    funValue[7] =0.125*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
-   funValue[8] =0.25*(1.-gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]);
-   funValue[9] =0.25*(1.-gc[0]*gc[0])*(1.+gc[1])*(1.-gc[2]);
-   funValue[10]=0.25*(1.+gc[0])*(1.-gc[1]*gc[1])*(1.-gc[2]);
-   funValue[11]=0.25*(1.-gc[0]*gc[0])*(1.-gc[1])*(1.-gc[2]);
-   funValue[12]=0.25*(1.-gc[0])*(1.-gc[1]*gc[1])*(1.+gc[2]);
-   funValue[13]=0.25*(1.-gc[0]*gc[0])*(1.+gc[1])*(1.+gc[2]);
-   funValue[14]=0.25*(1.+gc[0])*(1.-gc[1]*gc[1])*(1.+gc[2]);
-   funValue[15]=0.25*(1.-gc[0]*gc[0])*(1.-gc[1])*(1.+gc[2]);
-   funValue[16]=0.25*(1.-gc[0])*(1.-gc[1])*(1.-gc[2]*gc[2]);
-   funValue[17]=0.25*(1.-gc[0])*(1.+gc[1])*(1.-gc[2]*gc[2]);
-   funValue[18]=0.25*(1.+gc[0])*(1.+gc[1])*(1.-gc[2]*gc[2]);
-   funValue[19]=0.25*(1.+gc[0])*(1.-gc[1])*(1.-gc[2]*gc[2]);
+   funValue[8] =0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
+   funValue[9] =0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]-1.);
+   funValue[10]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
+   funValue[11]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]-1.);
+   funValue[12]=0.25*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
+   funValue[13]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*gc[2]*(gc[2]+1.);
+   funValue[14]=0.25*gc[0]*(gc[0]+1.)*(1.-gc[1]*gc[1])*gc[2]*(gc[2]+1.);
+   funValue[15]=0.25*(1.-gc[0]*gc[0])*gc[1]*(gc[1]-1.)*gc[2]*(gc[2]+1.);
+   funValue[16]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
+   funValue[17]=0.25*gc[0]*(gc[0]-1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
+   funValue[18]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
+   funValue[19]=0.25*gc[0]*(gc[0]+1.)*gc[1]*(gc[1]-1.)*(1.-gc[2]*gc[2]);
    funValue[20]=0.5*(1.-gc[0]*gc[0])*(1.-gc[1]*gc[1])*gc[2]*(gc[2]-1.);
    funValue[21]=0.5*gc[0]*(gc[0]-1.)*(1.-gc[1]*gc[1])*(1.-gc[2]*gc[2]);
    funValue[22]=0.5*(1.-gc[0]*gc[0])*gc[1]*(gc[1]+1.)*(1.-gc[2]*gc[2]);
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index 026fb199e..6bc2ddd4f 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -14439,6 +14439,33 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertTrue(d[0:].isEqual(DataArrayDouble([]),1e-12))
         pass
 
+    def testSwig2Hexa27GP1(self):
+        """ This test focused on shape functions of hexa27.
+        """
+        coo=DataArrayDouble([[0.,2.,2.],[0.,0.,2.],[2.,0.,2.],[2.,2.,2.],[0.,2.,0.],[0.,0.,0.],[2.,0.,0.],[2.,2.,0.], [0.,1.,2.],[1.,0.,2.],[2.,1.,2.],[1.,2.,2.], [0.,1.,0.],[1.,0.,0.],[2.,1.,0.],[1.,2.,0.], [0.,2.,1.],[0.,0.,1.],[2.,0.,1.],[2.,2.,1.], [1.,1.,2.], [0.,1.,1.],[1.,0.,1.],[2.,1.,1.],[1.,2.,1.], [1.,1.,0.], [1.,1.,1.]])
+        m=MEDCouplingUMesh("mesh",3) ; m.setCoords(coo)
+        m.allocateCells()
+        # the cell description is exactly those described in the description of HEXA27 in MED file 3.0.7 documentation
+        m.insertNextCell(NORM_HEXA27,[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26])
+        refCoo=[-1.,-1.,-1.,-1.,1.,-1.,1.,1.,-1.,1.,-1.,-1.,-1.,-1.,1.,-1.,1.,1.,1.,1.,1.,1.,-1.,1.,-1.,0.,-1.,0.,1.,-1.,1.,0.,-1.,0.,-1.,-1.,-1.,0.,1.,0.,1.,1.,1.,0.,1.,0.,-1.,1.,-1.,-1.,0.,-1.,1.,0.,1.,1.,0.,1.,-1.,0.,0.,0.,-1.,-1.,0.,0.,0.,1.,0.,1.,0.,0.,0.,-1.,0.,0.,0.,1.,0.,0.,0.]
+        weights=[0.1714677640603571,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.43895747599451346,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.43895747599451346,0.27434842249657115,0.43895747599451346,0.7023319615912209,0.43895747599451346,0.27434842249657115,0.43895747599451346,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.43895747599451346,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.1714677640603571]
+        gCoords=[-0.774596669241483,-0.774596669241483,-0.774596669241483,-0.774596669241483,-0.774596669241483,0.0,-0.774596669241483,-0.774596669241483,0.774596669241483,-0.774596669241483,0.0,-0.774596669241483,-0.774596669241483,0.0,0.0,-0.774596669241483,0.0,0.774596669241483,-0.774596669241483,0.774596669241483,-0.774596669241483,-0.774596669241483,0.774596669241483,0.0,-0.774596669241483,0.774596669241483,0.774596669241483,0.0,-0.774596669241483,-0.774596669241483,0.0,-0.774596669241483,0.0,0.0,-0.774596669241483,0.774596669241483,0.0,0.0,-0.774596669241483,0.0,0.0,0.0,0.0,0.0,0.774596669241483,0.0,0.774596669241483,-0.774596669241483,0.0,0.774596669241483,0.0,0.0,0.774596669241483,0.774596669241483,0.774596669241483,-0.774596669241483,-0.774596669241483,0.774596669241483,-0.774596669241483,0.0,0.774596669241483,-0.774596669241483,0.774596669241483,0.774596669241483,0.0,-0.774596669241483,0.774596669241483,0.0,0.0,0.774596669241483,0.0,0.774596669241483,0.774596669241483,0.774596669241483,-0.774596669241483,0.774596669241483,0.774596669241483,0.0,0.774596669241483,0.774596669241483,0.774596669241483]
+        fGauss=MEDCouplingFieldDouble(ON_GAUSS_PT) ; fGauss.setName("fGauss")
+        fGauss.setMesh(m)
+        fGauss.setGaussLocalizationOnType(NORM_HEXA27,refCoo,gCoords,weights)
+        arr=DataArrayDouble(fGauss.getNumberOfTuplesExpected()) ; arr.iota()
+        fGauss.setArray(arr)
+        arrOfDisc=fGauss.getLocalizationOfDiscr()
+        # the test is here
+        self.assertTrue(arrOfDisc.isEqual(DataArrayDouble([0.2254033307585172,1.7745966692414836,1.7745966692414834,0.22540333075851715,1.7745966692414834,1.,0.22540333075851715,1.7745966692414836,0.22540333075851715,0.22540333075851715,1.,1.7745966692414834,0.2254033307585171,1.,1.,0.22540333075851715,1.0000000000000002,0.2254033307585171,0.22540333075851715,0.22540333075851715,1.7745966692414838,0.22540333075851715,0.22540333075851715,1.,0.22540333075851715,0.22540333075851715,0.22540333075851715,1.,1.7745966692414832,1.7745966692414834,1.,1.774596669241483,1.,1.0000000000000002,1.7745966692414832,0.22540333075851712,1.,1.,1.774596669241483,1.,1.,1.,1.,1.,0.2254033307585171,1.,0.22540333075851715,1.7745966692414834,1.,0.2254033307585171,1.,1.0000000000000002,0.22540333075851715,0.2254033307585171,1.7745966692414834,1.7745966692414834,1.7745966692414836,1.7745966692414832,1.7745966692414834,1.0000000000000002,1.7745966692414834,1.7745966692414836,0.22540333075851712,1.7745966692414832,1.,1.7745966692414834,1.774596669241483,1.,1.,1.7745966692414832,1.0000000000000002,0.22540333075851712,1.7745966692414836,0.22540333075851715,1.7745966692414836,1.7745966692414832,0.22540333075851715,1.,1.7745966692414836,0.22540333075851715,0.22540333075851715],27,3),1e-12))
+        #
+        weights=27*[1]
+        gCoords=refCoo
+        fGauss.setGaussLocalizationOnType(NORM_HEXA27,refCoo,gCoords,weights)
+        arrOfDisc2=fGauss.getLocalizationOfDiscr()
+        self.assertTrue(arrOfDisc2.isEqual(coo,1e-12))
+        pass
+
     def setUp(self):
         pass
     pass
diff --git a/src/MEDLoader/Swig/MEDLoaderTest4.py b/src/MEDLoader/Swig/MEDLoaderTest4.py
index b176a7c7b..5a39364d2 100644
--- a/src/MEDLoader/Swig/MEDLoaderTest4.py
+++ b/src/MEDLoader/Swig/MEDLoaderTest4.py
@@ -4150,6 +4150,85 @@ class MEDLoaderTest4(unittest.TestCase):
         vExp=DataArrayDouble([-1.1,-3.1,5,4,3,2,1,0])
         self.assertTrue(v.isEqual(vExp,1e-12))
         pass
+
+    def test29(self):
+        """ This test focused on HEXA27 cell for which the MED numbering is not equal to the VTK numbering. So here the HEXA27 cell is those in MED file documentation (reference element).
+        """
+        fname="ForMEDReader29.med"
+        coo=DataArrayDouble([[0.,2.,2.],[0.,0.,2.],[2.,0.,2.],[2.,2.,2.],[0.,2.,0.],[0.,0.,0.],[2.,0.,0.],[2.,2.,0.], [0.,1.,2.],[1.,0.,2.],[2.,1.,2.],[1.,2.,2.], [0.,1.,0.],[1.,0.,0.],[2.,1.,0.],[1.,2.,0.], [0.,2.,1.],[0.,0.,1.],[2.,0.,1.],[2.,2.,1.], [1.,1.,2.], [0.,1.,1.],[1.,0.,1.],[2.,1.,1.],[1.,2.,1.], [1.,1.,0.], [1.,1.,1.]])
+        m=MEDCouplingUMesh("mesh",3) ; m.setCoords(coo)
+        m.allocateCells()
+        # MED = [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]
+        # VTK = [0,1,2,3,4,5,6,7, 8,9,10,11,12,13,14,15,16,17,18,19,24,22,21,23,20,25,26]
+        m.insertNextCell(NORM_HEXA27,[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26])
+        fCell=MEDCouplingFieldDouble(ON_CELLS) ; fCell.setName("fCell")
+        arrCell=DataArrayDouble([7.]) ; arrCell.setInfoOnComponent(0,"smth") ; fCell.setArray(arrCell)
+        fCell.setMesh(m)
+        MEDLoader.WriteField(fname,fCell,True)
+        refCoo=[-1.,-1.,-1.,-1.,1.,-1.,1.,1.,-1.,1.,-1.,-1.,-1.,-1.,1.,-1.,1.,1.,1.,1.,1.,1.,-1.,1.,-1.,0.,-1.,0.,1.,-1.,1.,0.,-1.,0.,-1.,-1.,-1.,0.,1.,0.,1.,1.,1.,0.,1.,0.,-1.,1.,-1.,-1.,0.,-1.,1.,0.,1.,1.,0.,1.,-1.,0.,0.,0.,-1.,-1.,0.,0.,0.,1.,0.,1.,0.,0.,0.,-1.,0.,0.,0.,1.,0.,0.,0.]
+        weights=[0.1714677640603571,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.43895747599451346,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.43895747599451346,0.27434842249657115,0.43895747599451346,0.7023319615912209,0.43895747599451346,0.27434842249657115,0.43895747599451346,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.43895747599451346,0.27434842249657115,0.1714677640603571,0.27434842249657115,0.1714677640603571]
+        gCoords=[-0.774596669241483,-0.774596669241483,-0.774596669241483,-0.774596669241483,-0.774596669241483,0.0,-0.774596669241483,-0.774596669241483,0.774596669241483,-0.774596669241483,0.0,-0.774596669241483,-0.774596669241483,0.0,0.0,-0.774596669241483,0.0,0.774596669241483,-0.774596669241483,0.774596669241483,-0.774596669241483,-0.774596669241483,0.774596669241483,0.0,-0.774596669241483,0.774596669241483,0.774596669241483,0.0,-0.774596669241483,-0.774596669241483,0.0,-0.774596669241483,0.0,0.0,-0.774596669241483,0.774596669241483,0.0,0.0,-0.774596669241483,0.0,0.0,0.0,0.0,0.0,0.774596669241483,0.0,0.774596669241483,-0.774596669241483,0.0,0.774596669241483,0.0,0.0,0.774596669241483,0.774596669241483,0.774596669241483,-0.774596669241483,-0.774596669241483,0.774596669241483,-0.774596669241483,0.0,0.774596669241483,-0.774596669241483,0.774596669241483,0.774596669241483,0.0,-0.774596669241483,0.774596669241483,0.0,0.0,0.774596669241483,0.0,0.774596669241483,0.774596669241483,0.774596669241483,-0.774596669241483,0.774596669241483,0.774596669241483,0.0,0.774596669241483,0.774596669241483,0.774596669241483]
+        fGauss=MEDCouplingFieldDouble(ON_GAUSS_PT) ; fGauss.setName("fGauss")
+        fGauss.setMesh(m)
+        fGauss.setGaussLocalizationOnType(NORM_HEXA27,refCoo,gCoords,weights)
+        arrGauss=DataArrayDouble(fGauss.getNumberOfTuplesExpected()) ; arrGauss.setInfoOnComponent(0,"gaussc") ; arrGauss.iota()
+        fGauss.setArray(arrGauss)
+        MEDLoader.WriteFieldUsingAlreadyWrittenMesh(fname,fGauss)
+        ########## GO for reading in MEDReader,by not loading all. Mesh is fully loaded but not fields values
+        ms=MEDFileMeshes(fname)
+        fields=MEDFileFields(fname,False)
+        fields.removeFieldsWithoutAnyTimeStep()
+        fields_per_mesh=[fields.partOfThisLyingOnSpecifiedMeshName(meshName) for meshName in ms.getMeshesNames()]
+        allFMTSLeavesToDisplay=[]
+        for fields in fields_per_mesh:
+            allFMTSLeavesToDisplay2=[]
+            for fmts in fields:
+                allFMTSLeavesToDisplay2+=fmts.splitDiscretizations()
+                pass
+            allFMTSLeavesToDisplay.append(allFMTSLeavesToDisplay2)
+            pass
+        self.assertEqual(len(allFMTSLeavesToDisplay),1)
+        self.assertEqual(len(allFMTSLeavesToDisplay[0]),2)
+        allFMTSLeavesPerTimeSeries=MEDFileAnyTypeFieldMultiTS.SplitIntoCommonTimeSeries(sum(allFMTSLeavesToDisplay,[]))
+        self.assertEqual(len(allFMTSLeavesPerTimeSeries),1)
+        self.assertEqual(len(allFMTSLeavesPerTimeSeries[0]),2)
+        allFMTSLeavesPerCommonSupport1=MEDFileAnyTypeFieldMultiTS.SplitPerCommonSupport(allFMTSLeavesToDisplay[0],ms[ms.getMeshesNames()[0]])
+        self.assertEqual(len(allFMTSLeavesPerCommonSupport1),1)
+        self.assertEqual(len(allFMTSLeavesPerCommonSupport1[0][0]),2)
+        #
+        mst=MEDFileMeshStruct.New(ms[0])
+        #
+        fcscp=allFMTSLeavesPerCommonSupport1[0][1]
+        mml=fcscp.buildFromScratchDataSetSupport(0,fields)
+        mml2=mml.prepare()
+        self.assertTrue(isinstance(mml2,MEDUMeshMultiLev))
+        ncc,a0,a1,a2,a3,a4,a5=mml2.buildVTUArrays()
+        self.assertTrue(a0.isEqual(coo,1e-12))
+        self.assertTrue(a1.isEqual(DataArrayByte([29])))
+        self.assertTrue(a2.isEqual(DataArrayInt([0])))
+        # the connectivity must be not a iota as declared in m.insertNextCell
+        self.assertTrue(a3.isEqual(DataArrayInt([27,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,24,22,21,23,20,25,26])))# the test is on this line to check that connectivity has been processed for HEXA27
+        self.assertTrue(a4 is None)
+        self.assertTrue(a5 is None)
+        ffCell=allFMTSLeavesPerCommonSupport1[0][0][0][0]
+        fsst=MEDFileField1TSStructItem.BuildItemFrom(ffCell,mst)
+        ffCell.loadArraysIfNecessary()
+        v=mml2.buildDataArray(fsst,fields,ffCell.getUndergroundDataArray())
+        self.assertEqual(v.getHiddenCppPointer(),ffCell.getUndergroundDataArray().getHiddenCppPointer())
+        self.assertEqual(ffCell.getName(),"fCell")
+        self.assertTrue(v.isEqual(arrCell,1e-12)) ; self.assertTrue(v.isEqualWithoutConsideringStr(DataArrayDouble([7.]),1e-12)) ; self.assertEqual(v.getInfoOnComponents(),["smth"])
+        del ffCell
+        #
+        ffGauss=allFMTSLeavesPerCommonSupport1[0][0][1][0]
+        fsst=MEDFileField1TSStructItem.BuildItemFrom(ffGauss,mst)
+        ffGauss.loadArraysIfNecessary()
+        v=mml2.buildDataArray(fsst,fields,ffGauss.getUndergroundDataArray())
+        self.assertEqual(v.getHiddenCppPointer(),ffGauss.getUndergroundDataArray().getHiddenCppPointer())
+        self.assertEqual(ffGauss.getName(),"fGauss")
+        self.assertTrue(v.isEqual(arrGauss,1e-12)) ; self.assertTrue(v.isEqualWithoutConsideringStr(DataArrayDouble(range(27)),1e-12)) ; self.assertEqual(v.getInfoOnComponents(),["gaussc"])
+        ffGauss=allFMTSLeavesPerCommonSupport1[0][0][1][0]
+        pass
+
     pass
 
 unittest.main()
-- 
2.39.2