From: Anthony Geay Date: Tue, 31 Jan 2023 08:27:55 +0000 (+0100) Subject: [EDF26877] : management of computation of measure field on field on Gauss Point in... X-Git-Tag: V9_11_0a1~14 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=9abc40a840f69fd7f07b356cd31876b64dc81e14;p=tools%2Fmedcoupling.git [EDF26877] : management of computation of measure field on field on Gauss Point in the special case where SpaceDim is not equal to MeshDim --- diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx index 05435efe3..72c49b896 100644 --- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx +++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx @@ -848,6 +848,14 @@ void GaussInfo::seg2aInit() funValue[0] = 0.5*(1.0 - gc[0]); funValue[1] = 0.5*(1.0 + gc[0]); SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + + devFunValue[0] = -0.5; + + devFunValue[1] = 0.5; + + DEV_SHAPE_FUN_MACRO_END; } void GaussInfo::seg2bInit() @@ -865,6 +873,14 @@ void GaussInfo::seg2bInit() funValue[0] = 1.0 - gc[0]; funValue[1] = gc[0]; SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + + devFunValue[0] = -1.0 ; + + devFunValue[1] = 1.0 ; + + DEV_SHAPE_FUN_MACRO_END; } /*! diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx index 63080077f..f363e2122 100644 --- a/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx +++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.hxx @@ -48,6 +48,7 @@ namespace INTERP_KERNEL void assign(mcIdType newn, mcIdType newm, const T &a); ~DenseMatrixT(); T determinant() const; + T toJacobian() const; }; using DenseMatrix = DenseMatrixT; diff --git a/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx index 2a262c31b..3ae1aa2c3 100644 --- a/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx +++ b/src/INTERP_KERNEL/InterpKernelDenseMatrix.txx @@ -21,6 +21,9 @@ #include "InterpKernelDenseMatrix.hxx" #include "InterpKernelException.hxx" +#include "VectorUtils.hxx" + +#include namespace INTERP_KERNEL { @@ -151,10 +154,33 @@ namespace INTERP_KERNEL template T DenseMatrixT::determinant() const { + if(nn==1 && mm==1) + return v[0][0]; if(nn==2 && mm==2) return Determinant22(v[0]); if(nn==3 && mm==3) return Determinant33(v[0]); - THROW_IK_EXCEPTION("DenseMatrixT::determinant : only 2x2 and 3x3 implemented !"); + THROW_IK_EXCEPTION("DenseMatrixT::determinant : only 1x1, 2x2 and 3x3 implemented !"); + } + + template + T DenseMatrixT::toJacobian() const + { + if(nn == mm) + return determinant(); + const T *vPtr(this->v[0]); + if(nn==3 && mm==1) + return norm(vPtr); + if(nn==2 && mm==1) + return std::sqrt(vPtr[0]*vPtr[0] + vPtr[1]*vPtr[1]); + if(nn==3 && mm==2) + { + T tmp[3]; + T VA[3] = {v[0][0],v[1][0],v[2][0]}; + T VB[3] = {v[0][1],v[1][1],v[2][1]}; + cross(VA,VB,tmp); + return norm(tmp); + } + THROW_IK_EXCEPTION("DenseMatrixT::toJacobian : only 1x1, 2x1, 3x1, 3x2, 2x2 and 3x3 implemented !"); } } diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index ebf44b23d..d35dccf51 100755 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -1715,8 +1715,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con MCAuto umesh(mesh->buildUnstructured()); const double *coordsOfMesh( umesh->getCoords()->begin() ); auto spaceDim(mesh->getSpaceDimension()); - if( spaceDim != mesh->getMeshDimension()) - THROW_IK_EXCEPTION("MEDCouplingFieldDiscretizationGauss::getMeasureField : only implemented when space dimension and mesh dimension are equal !"); + auto meshDim(mesh->getMeshDimension()); MCAuto ret=MEDCouplingFieldDouble::New(ON_GAUSS_PT); ret->setMesh(mesh); ret->setDiscretization(const_cast(this)); @@ -1752,19 +1751,19 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con std::vector ptsInCell; ptsInCell.reserve(conn.size()*loc.getDimension()); std::for_each( conn.cbegin(), conn.cend(), [spaceDim,coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*spaceDim,coordsOfMesh+(c+1)*spaceDim); } ); std::size_t nbPtsInCell(ptsInCell.size()/spaceDim); - INTERP_KERNEL::DenseMatrix jacobian(spaceDim,spaceDim); + INTERP_KERNEL::DenseMatrix jacobian(spaceDim,meshDim); MCAuto shapeFunc = loc.getDerivativeOfShapeFunctionValues(); for(mcIdType iGPt = 0 ; iGPt < nbOfGaussPt ; ++iGPt) { for(auto i = 0 ; i < spaceDim ; ++i) - for(auto j = 0 ; j < spaceDim ; ++j) + for(auto j = 0 ; j < meshDim ; ++j) { double res = 0.0; for( std::size_t k = 0 ; k < nbPtsInCell ; ++k ) - res += ptsInCell[spaceDim*k+i] * shapeFunc->getIJ(iGPt,spaceDim*k+j); + res += ptsInCell[spaceDim*k+i] * shapeFunc->getIJ(iGPt,meshDim*k+j); jacobian[ i ][ j ] = res; } - arrPtr[offsetPtr[*cellId]+iGPt]=std::abs( jacobian.determinant() )*loc.getWeight(FromIdType(iGPt)); + arrPtr[offsetPtr[*cellId]+iGPt]=std::abs( jacobian.toJacobian() )*loc.getWeight(FromIdType(iGPt)); } } } diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index d8815278c..89c249da6 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1175,6 +1175,36 @@ class MEDCouplingBasicsTest7(unittest.TestCase): self.assertTrue( ret.isEqual(DataArrayInt([1, 1, 0, 1, 2, 1, 0, 1, 2]) ) ) pass + def testMeasureOnGaussPtMeshDimNotEqualSpaceDim0(self): + """ + [EDF26877] : This test focuses on computation of measure field on field on Gauss Point in the special case where SpaceDim + are not eqaul to the meshDim. + """ + seg2 = MEDCouplingUMesh("mesh",1) + seg2.setCoords(DataArrayDouble([(3,3),(4,4)])) + seg2.allocateCells() + seg2.insertNextCell(NORM_SEG2,[0,1]) + fff=MEDCouplingFieldDouble.New(ON_GAUSS_PT) ; fff.setName("CH1RB") ; fff.setNature(IntensiveMaximum) + fff.setMesh(seg2) + fff.setGaussLocalizationOnCells([0], [0.,1.], [0.333333333333333], [1.0]) + disc = fff.getDiscretization() + # spaceDim = 2 meshDim = 1 + self.assertTrue( disc.getMeasureField(seg2,True).getArray().isEqual(DataArrayDouble([sqrt(2.0)]),1e-10) ) + # spaceDim = 3 meshDim = 1 + seg2.setCoords(DataArrayDouble([(3,3,3),(4,4,4)])) + self.assertTrue( disc.getMeasureField(seg2,True).getArray().isEqual(DataArrayDouble([sqrt(3.0)]),1e-10) ) + # spaceDim = 3 meshDim = 2 + tri = MEDCouplingUMesh("mesh",2) + tri.setCoords( DataArrayDouble([(0,0,0),(1,1,0),(2,2,2)]) ) + tri.allocateCells() + tri.insertNextCell(NORM_TRI3,[0,1,2]) + fff=MEDCouplingFieldDouble.New(ON_GAUSS_PT) ; fff.setName("CH1RB") ; fff.setNature(IntensiveMaximum) + fff.setMesh(tri) + fff.setGaussLocalizationOnCells(list(range(0, 1)), [0., 0., 1., 0., 0., 1.], [0.3333333333333333, 0.3333333333333333], [0.5]) + disc = fff.getDiscretization() + self.assertTrue( disc.getMeasureField(tri,True).getArray().isEqual( tri.getMeasureField(True).getArray(), 1e-10) ) + pass + pass if __name__ == '__main__':