From cca438f59fd43e1a9f0d980bcac870f95bad81a5 Mon Sep 17 00:00:00 2001 From: Anthony Geay Date: Fri, 4 Aug 2023 09:39:13 +0200 Subject: [PATCH] [EDF28170] : management of structure elements pipes using SEG4 elements. Management of Gauss Points on SEG4 cell types --- .../GaussPoints/InterpKernelGaussCoords.cxx | 46 +++++++++ .../GaussPoints/InterpKernelGaussCoords.hxx | 2 + .../MEDCouplingBasicsTest7.py | 2 +- src/MEDLoader/MEDFileBlowStrEltUp.cxx | 99 ++++++++++++++++--- 4 files changed, 137 insertions(+), 12 deletions(-) diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx index 85752d6ba..12cb24902 100644 --- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx +++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx @@ -35,6 +35,8 @@ const double GaussInfo::SEG2B_REF[2]={0., 1.0}; const double GaussInfo::SEG3_REF[3]={-1.0, 1.0, 0.0}; +const double GaussInfo::SEG4_REF[4]={-1.0, 1.0, -0.3333333333333333, 0.3333333333333333}; + const double GaussInfo::TRIA3A_REF[6]={-1.0, 1.0, -1.0, -1.0, 1.0, -1.0}; const double GaussInfo::TRIA3B_REF[6]={0.0, 0.0, 1.0, 0.0, 0.0, 1.0}; @@ -471,6 +473,8 @@ std::vector GaussInfo::GetDefaultReferenceCoordinatesOf(NormalizedCellTy return std::vector(SEG2A_REF,SEG2A_REF+sizeof(SEG2A_REF)/sizeof(double)); case INTERP_KERNEL::NORM_SEG3: return std::vector(SEG3_REF,SEG3_REF+sizeof(SEG3_REF)/sizeof(double)); + case INTERP_KERNEL::NORM_SEG4: + return std::vector(SEG4_REF,SEG4_REF+sizeof(SEG4_REF)/sizeof(double)); case INTERP_KERNEL::NORM_TRI3: return std::vector(TRIA3A_REF,TRIA3A_REF+sizeof(TRIA3A_REF)/sizeof(double)); case INTERP_KERNEL::NORM_TRI6: @@ -578,6 +582,14 @@ void GaussInfo::initLocalInfo() CHECK_MACRO; break; + case NORM_SEG4: + _my_local_ref_dim = 1; + _my_local_nb_ref = 4; + seg4Init(); + aSatify = isSatisfy(); + CHECK_MACRO; + break; + case NORM_TRI3: _my_local_ref_dim = 2; _my_local_nb_ref = 3; @@ -912,6 +924,40 @@ void GaussInfo::seg3Init() DEV_SHAPE_FUN_MACRO_END; } +/*! + * Init Segment 4 Reference coordinates ans Shape function. + */ +void GaussInfo::seg4Init() +{ + LOCAL_COORD_MACRO_BEGIN; + case 0: + coords[0] = SEG4_REF[0]; + break; + case 1: + coords[0] = SEG4_REF[1]; + break; + case 2: + coords[0] = SEG4_REF[2]; + break; + case 3: + coords[0] = SEG4_REF[3]; + LOCAL_COORD_MACRO_END; + + SHAPE_FUN_MACRO_BEGIN; + funValue[0] = 9.0/16.0 * (1-gc[0])*(gc[0]+1.0/3.0)*(gc[0]-1.0/3.0); + funValue[1] = -9.0/16.0 * (1+gc[0])*(1.0/3.0-gc[0])*(gc[0]+1.0/3.0); + funValue[2] = 27.0/16.0 * (gc[0]-1)*(gc[0]+1)*(gc[0]-1.0/3.0); + funValue[3] = -27.0/16.0 * (gc[0]-1)*(gc[0]+1)*(gc[0]+1.0/3.0); + SHAPE_FUN_MACRO_END; + + DEV_SHAPE_FUN_MACRO_BEGIN; + devFunValue[0] = 9.0/16.0 * (-3.0*gc[0]*gc[0]+2.0*gc[0]+1.0/9.0); + devFunValue[1] = 9.0/16.0 * (3.0*gc[0]*gc[0]+2.0*gc[0]-1.0/9.0); + devFunValue[2] = 27.0/16.0 * (3.0*gc[0]*gc[0]-2.0/3.0*gc[0]-1.0); + devFunValue[3] = -27.0/16.0 * (3.0*gc[0]*gc[0]+2.0/3.0*gc[0]-1.0); + DEV_SHAPE_FUN_MACRO_END; +} + /*! * Init Triangle Reference coordinates ans Shape function. * Case A. diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx index a1c21fd59..948877a5e 100644 --- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx +++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx @@ -69,6 +69,7 @@ namespace INTERP_KERNEL static const double SEG2A_REF[2]; static const double SEG2B_REF[2]; static const double SEG3_REF[3]; + static const double SEG4_REF[4]; static const double TRIA3A_REF[6]; static const double TRIA3B_REF[6]; static const double TRIA6A_REF[12]; @@ -108,6 +109,7 @@ namespace INTERP_KERNEL void seg2aInit(); void seg2bInit(); void seg3Init(); + void seg4Init(); //2D void tria3aInit(); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py index c32bd34bd..4f3a298f1 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest7.py @@ -1174,7 +1174,7 @@ class MEDCouplingBasicsTest7(unittest.TestCase): # 1D cells vec = [0.64] - for gt in [NORM_SEG2,NORM_SEG3]: + for gt in [NORM_SEG2,NORM_SEG3,NORM_SEG4]: ref_coord = [list(elt) for elt in MEDCouplingGaussLocalization.GetDefaultReferenceCoordinatesOf(gt).getValuesAsTuple()] der_computed = GetDerivative(ref_coord,vec) diff --git a/src/MEDLoader/MEDFileBlowStrEltUp.cxx b/src/MEDLoader/MEDFileBlowStrEltUp.cxx index 1a0465ab5..1292c0c8d 100644 --- a/src/MEDLoader/MEDFileBlowStrEltUp.cxx +++ b/src/MEDLoader/MEDFileBlowStrEltUp.cxx @@ -340,6 +340,7 @@ private: static MCAuto BuildMeshFromAngleVrille(INTERP_KERNEL::NormalizedCellType gt, const DataArrayDouble *angleDeVrille, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs); static MCAuto BuildMeshFromEpaisseur(INTERP_KERNEL::NormalizedCellType gt, const DataArrayDouble *thickness, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs); static MCAuto BuildMeshPipeSEG3(const DataArrayDouble *angle, const DataArrayDouble *scale, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs); + static MCAuto BuildMeshPipeSEG4(const DataArrayDouble *angle, const DataArrayDouble *scale, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs); static MCAuto BuildMeshCommon(INTERP_KERNEL::NormalizedCellType gt, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs, MCAuto& ptsForLoc); static MCAuto BuildMeshFromStructure(INTERP_KERNEL::NormalizedCellType gt, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs); public: @@ -579,6 +580,78 @@ MCAuto LocInfo::BuildMeshPipeSEG3(const DataArrayDouble *angle, return resu; } +MCAuto LocInfo::BuildMeshPipeSEG4(const DataArrayDouble *angle, const DataArrayDouble *scale, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs) +{ + static const char MSG1[]="BuildMeshPipeSEG4 : not recognized pattern ! Send mail to anthony.geay@edf.fr with corresponding MED file !"; + MCAuto ptsForLoc; + MCAuto geoMesh(BuildMeshCommon(INTERP_KERNEL::NORM_SEG4,pfl,loc,zeStr,mesh,section,globs,ptsForLoc)); + mcIdType nbSecPts(section->getNumberOfNodes()),nbCells(geoMesh->getNumberOfCells()),nbg(ToIdType(loc.getGaussWeights().size())); + MCConstAuto zeAngle,zeScale; + if(!pfl.empty()) + { + const DataArrayIdType *pflArr(globs->getProfile(pfl)); + zeAngle=angle->selectByTupleIdSafe(pflArr->begin(),pflArr->end()); + zeScale=scale->selectByTupleIdSafe(pflArr->begin(),pflArr->end()); + } + else + { + zeAngle.takeRef(angle); + zeScale.takeRef(scale); + } + if(zeAngle->getNumberOfComponents()!=4 || zeScale->getNumberOfComponents()!=2 || nbg!=3) + throw INTERP_KERNEL::Exception(MSG1); + MCAuto dir; + { + MCAuto geoMesh2(geoMesh->deepCopy()); + geoMesh2->convertQuadraticCellsToLinear(); + dir=geoMesh2->buildDirectionVectorField(); + } + MCAuto rot(dir->getArray()->fromCartToSpher()); + std::size_t nbCompo(ptsForLoc->getNumberOfComponents()); + MCAuto secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.)); + { + const int TAB[3]={2,0,1}; + std::vector v(TAB,TAB+3); + secPts=secPts->keepSelectedComponents(v); + } + const double CENTER[3]={0.,0.,0.},AX0[3]={0.,0.,1.}; + double AX1[3]; AX1[2]=0.; + std::vector< MCAuto > arrs(nbCells*nbg); + for(mcIdType j=0;j p(secPts->deepCopy()); + double ang0(rot->getIJ(j,2)); + double rmin(zeScale->getIJ(j,0)),rmax(zeScale->getIJ(j,1)); + { + auto pt(p->rwBegin()); + for(int i=0;ibegin(),p->getPointer()); + AX1[0]=-sin(ang0); AX1[1]=cos(ang0);// rot Oy around OZ + double ang1(M_PI/2.-rot->getIJ(j,1)); + DataArrayDouble::Rotate3DAlg(CENTER,AX1,-ang1,nbSecPts,p->begin(),p->getPointer()); + for(int l=0;l<3;l++) + { + MCAuto p3(p->deepCopy()); + DataArrayDouble::Rotate3DAlg(CENTER,dir->getArray()->begin()+j*3,zeAngle->getIJ(j,l),nbSecPts,p3->begin(),p3->getPointer()); + MCAuto p2(p3->deepCopy()); + for(std::size_t k=0;kapplyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k); + arrs[j*nbg+l]=p2; + } + } + std::vector arrs2(VecAutoToVecOfCstPt(arrs)); + MCAuto resu(DataArrayDouble::Aggregate(arrs2)); + return resu; +} + MCAuto LocInfo::BuildMeshFromStructure(INTERP_KERNEL::NormalizedCellType gt, const std::string& pfl, const MEDFileFieldLoc& loc, const MEDFileEltStruct4Mesh *zeStr, const MEDFileUMesh *mesh, const MEDFileUMesh *section, const MEDFileFieldGlobsReal *globs) { static const char MSG1[]="BuildMeshFromStructure : not recognized pattern ! Send mail to anthony.geay@edf.fr with corresponding MED file !"; @@ -602,23 +675,27 @@ MCAuto LocInfo::BuildMeshFromStructure(INTERP_KERNEL::Normalize if(zeArr0.isNull() || zeArr1.isNull()) throw INTERP_KERNEL::Exception(MSG1); MCAuto zeArr00(DynamicCastSafe(zeArr0)),zeArr11(DynamicCastSafe(zeArr1)); + MCAuto angle,scale; + if(zeArr00->getName()==ANGLE) + angle=zeArr00; + if(zeArr00->getName()==SCALE) + scale=zeArr00; + if(zeArr11->getName()==ANGLE) + angle=zeArr11; + if(zeArr11->getName()==SCALE) + scale=zeArr11; + if(angle.isNull() || scale.isNull()) + throw INTERP_KERNEL::Exception(MSG1); switch(gt) { case INTERP_KERNEL::NORM_SEG3: { - MCAuto angle,scale; - if(zeArr00->getName()==ANGLE) - angle=zeArr00; - if(zeArr00->getName()==SCALE) - scale=zeArr00; - if(zeArr11->getName()==ANGLE) - angle=zeArr11; - if(zeArr11->getName()==SCALE) - scale=zeArr11; - if(angle.isNull() || scale.isNull()) - throw INTERP_KERNEL::Exception(MSG1); return BuildMeshPipeSEG3(angle,scale,pfl,loc,zeStr,mesh,section,globs); } + case INTERP_KERNEL::NORM_SEG4: + { + return BuildMeshPipeSEG4(angle,scale,pfl,loc,zeStr,mesh,section,globs); + } default: throw INTERP_KERNEL::Exception(MSG1); } -- 2.39.2