[EDF28170] : management of structure elements pipes using SEG4 elements. Management...

diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.cxx
index 85752d6ba5d33cd1de27ea0cd793fa422619cc2c..12cb24902829b38296aaa404a506a7e72207ea39 100644 (file)
--- 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<double> GaussInfo::GetDefaultReferenceCoordinatesOf(NormalizedCellTy
       return std::vector<double>(SEG2A_REF,SEG2A_REF+sizeof(SEG2A_REF)/sizeof(double));
     case INTERP_KERNEL::NORM_SEG3:
       return std::vector<double>(SEG3_REF,SEG3_REF+sizeof(SEG3_REF)/sizeof(double));
+    case INTERP_KERNEL::NORM_SEG4:
+      return std::vector<double>(SEG4_REF,SEG4_REF+sizeof(SEG4_REF)/sizeof(double));
     case INTERP_KERNEL::NORM_TRI3:
       return std::vector<double>(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 a1c21fd59deb3ccdda1c2aa7931109d53f4f9b8a..948877a5ee01ea113275178bdb3c99bb1d02513f 100644 (file)
--- 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 c32bd34bd78c19d370a4e308666900cf80e58e90..4f3a298f102d529223f155ca9c8354c7663210d9 100644 (file)
--- 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 1a0465ab5325b60b2cf9e11167f87b4dc8da1844..1292c0c8d05e99add6cf6334cec3306272c63993 100644 (file)
--- a/src/MEDLoader/MEDFileBlowStrEltUp.cxx
+++ b/src/MEDLoader/MEDFileBlowStrEltUp.cxx
@@ -340,6 +340,7 @@ private:
   static MCAuto<DataArrayDouble> 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<DataArrayDouble> 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<DataArrayDouble> 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<DataArrayDouble> 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<MEDCouplingUMesh> 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<DataArrayDouble>& ptsForLoc);
   static MCAuto<DataArrayDouble> 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<DataArrayDouble> LocInfo::BuildMeshPipeSEG3(const DataArrayDouble *angle,
   return resu;
 }
 
+MCAuto<DataArrayDouble> 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<DataArrayDouble> ptsForLoc;
+  MCAuto<MEDCouplingUMesh> 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<DataArrayDouble> 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<MEDCouplingFieldDouble> dir;
+  {
+    MCAuto<MEDCouplingUMesh> geoMesh2(geoMesh->deepCopy());
+    geoMesh2->convertQuadraticCellsToLinear();
+    dir=geoMesh2->buildDirectionVectorField();
+  }
+  MCAuto<DataArrayDouble> rot(dir->getArray()->fromCartToSpher());
+  std::size_t nbCompo(ptsForLoc->getNumberOfComponents());
+  MCAuto<DataArrayDouble> secPts(section->getCoords()->changeNbOfComponents(nbCompo,0.));
+  {
+    const int TAB[3]={2,0,1};
+    std::vector<std::size_t> 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<DataArrayDouble> > arrs(nbCells*nbg);
+  for(mcIdType j=0;j<nbCells;j++)
+    {
+      constexpr int DIM=3;
+      MCAuto<DataArrayDouble> 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;i<nbSecPts;i++)
+          {
+            auto nrm(sqrt(std::accumulate(pt,pt+DIM,0.,[](double sum, double v) { return sum+v*v; } )));
+            auto sca((rmin+2.*(nrm-0.5)*(rmax-rmin))/nrm);
+            std::for_each(pt,pt+3,[sca](double& val) { val*=sca; } );
+            std::advance(pt,DIM);
+          }
+      }
+      DataArrayDouble::Rotate3DAlg(CENTER,AX0,ang0,nbSecPts,p->begin(),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<DataArrayDouble> p3(p->deepCopy());
+          DataArrayDouble::Rotate3DAlg(CENTER,dir->getArray()->begin()+j*3,zeAngle->getIJ(j,l),nbSecPts,p3->begin(),p3->getPointer());
+          MCAuto<DataArrayDouble> p2(p3->deepCopy());
+          for(std::size_t k=0;k<nbCompo;k++)
+            p2->applyLin(1.,ptsForLoc->getIJ(j*nbg+l,k),k);
+          arrs[j*nbg+l]=p2;
+        }
+    }
+  std::vector<const DataArrayDouble *> arrs2(VecAutoToVecOfCstPt(arrs));
+  MCAuto<DataArrayDouble> resu(DataArrayDouble::Aggregate(arrs2));
+  return resu;
+}
+
 MCAuto<DataArrayDouble> 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<DataArrayDouble> LocInfo::BuildMeshFromStructure(INTERP_KERNEL::Normalize
       if(zeArr0.isNull() || zeArr1.isNull())
         throw INTERP_KERNEL::Exception(MSG1);
       MCAuto<DataArrayDouble> zeArr00(DynamicCastSafe<DataArray,DataArrayDouble>(zeArr0)),zeArr11(DynamicCastSafe<DataArray,DataArrayDouble>(zeArr1));
+      MCAuto<DataArrayDouble> 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<DataArrayDouble> 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);
       }