New references checked for PYRA5 due to modification of shape function in fd8dbd34764...

[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingGaussLocalization.cxx

diff --git a/src/MEDCoupling/MEDCouplingGaussLocalization.cxx b/src/MEDCoupling/MEDCouplingGaussLocalization.cxx
index ebf9cc4a0a44a53167824d4b241fbcafe108f184..5f1920d96dcefdedc82dffcafddb06bec53eca27 100644 (file)
--- a/src/MEDCoupling/MEDCouplingGaussLocalization.cxx
+++ b/src/MEDCoupling/MEDCouplingGaussLocalization.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2020  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2022  CEA/DEN, EDF R&D
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -34,29 +34,18 @@ using namespace MEDCoupling;
 
 MEDCouplingGaussLocalization::MEDCouplingGaussLocalization(INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
                                                            const std::vector<double>& gsCoo, const std::vector<double>& w)
-try:_type(type),_ref_coord(refCoo),_gauss_coord(gsCoo),_weight(w)
+:_type(type),_ref_coord(refCoo),_gauss_coord(gsCoo),_weight(w)
 {
+  // Will potentially throw (and then release memory for above objects _ref_coord etc ...)
   checkConsistencyLight();
 }
-catch(INTERP_KERNEL::Exception& e)
-{
-    _type=INTERP_KERNEL::NORM_ERROR;
-    _ref_coord.clear();
-    _gauss_coord.clear();
-    _weight.clear();
-    throw e;
-}
 
 MEDCouplingGaussLocalization::MEDCouplingGaussLocalization(INTERP_KERNEL::NormalizedCellType typ)
-try:_type(typ)
+:_type(typ)
 {
+  // Will potentially throw
   INTERP_KERNEL::CellModel::GetCellModel(_type);
 }
-catch(INTERP_KERNEL::Exception& e)
-{
-    _type=INTERP_KERNEL::NORM_ERROR;
-    throw e;
-}
 
 void MEDCouplingGaussLocalization::setType(INTERP_KERNEL::NormalizedCellType typ)
 {
@@ -87,15 +76,20 @@ void MEDCouplingGaussLocalization::checkConsistencyLight() const
 int MEDCouplingGaussLocalization::getDimension() const
 {
   if(_weight.empty())
-    return -1;
+    THROW_IK_EXCEPTION("getDimension : weight is empty !");
   return (int)_gauss_coord.size()/(int)_weight.size();
 }
 
 int MEDCouplingGaussLocalization::getNumberOfPtsInRefCell() const
 {
-  int dim=getDimension();
-  if(dim==0)
-    return -1;
+  if(_gauss_coord.empty())
+  {
+    if(!_weight.empty())
+      THROW_IK_EXCEPTION("getNumberOfPtsInRefCell : gauss_coords are empty whereas weights are not empty !");
+    const INTERP_KERNEL::CellModel& cm = INTERP_KERNEL::CellModel::GetCellModel(_type);
+    return ((int)_ref_coord.size()) / ((int)cm.getDimension());
+  }
+  int dim( getDimension() );
   return (int)_ref_coord.size()/dim;
 }
 
@@ -149,7 +143,7 @@ double MEDCouplingGaussLocalization::getGaussCoord(int gaussPtIdInCell, int comp
   return _gauss_coord[gaussPtIdInCell*dim+comp];
 }
 
-double MEDCouplingGaussLocalization::getWeight(int gaussPtIdInCell, double newVal) const
+double MEDCouplingGaussLocalization::getWeight(int gaussPtIdInCell) const
 {
   checkCoherencyOfRequest(gaussPtIdInCell,0);
   return _weight[gaussPtIdInCell];
@@ -224,7 +218,6 @@ MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::localizePtsInRefCooForEach
   double *retPtr(ret->getPointer());
   if(dim!=ToIdType(ptsInRefCoo->getNumberOfComponents()))
     throw INTERP_KERNEL::Exception("MEDCouplingGaussLocalization::localizePtsInRefCooForEachCell : number of components of input coo is not equal to dim of element !");
-  const std::vector<double>& wg(getWeights());
   INTERP_KERNEL::GaussCoords calculator;
   calculator.addGaussInfo(typ,dim, ptsInRefCoo->begin(),nbPts,&_ref_coord[0],getNumberOfPtsInRefCell());
   //
@@ -253,6 +246,46 @@ MCAuto<MEDCouplingUMesh> MEDCouplingGaussLocalization::buildRefCell() const
   return MCAuto<MEDCouplingUMesh>(ret->buildUnstructured());
 }
 
+/*!
+ * This method returns an array containing for each Gauss Points in \a this, function values relative to the points of the
+ * reference cell. Number of components of returned array is equal to the number of points of the reference cell.
+ */
+MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::getShapeFunctionValues() const
+{
+  MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+  int nbGaussPt(getNumberOfGaussPt()),nbPtsRefCell(getNumberOfPtsInRefCell()),dim(getDimension());
+  ret->alloc(nbGaussPt,nbPtsRefCell);
+  double *retPtr(ret->getPointer());
+  for(int iGaussPt = 0 ; iGaussPt < nbGaussPt ; ++iGaussPt)
+  {
+    std::vector<double> curGaussPt(_gauss_coord.begin()+iGaussPt*dim,_gauss_coord.begin()+(iGaussPt+1)*dim);
+    INTERP_KERNEL::GaussInfo gi(_type,curGaussPt,1,_ref_coord,nbPtsRefCell);
+    gi.initLocalInfo();
+    const double *funcVal( gi.getFunctionValues(0) );
+    std::copy(funcVal,funcVal+nbPtsRefCell,retPtr);
+    retPtr += nbPtsRefCell;
+  }
+  return ret;
+}
+
+MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::getDerivativeOfShapeFunctionValues() const
+{
+  MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+  int nbGaussPt(getNumberOfGaussPt()),nbPtsRefCell(getNumberOfPtsInRefCell()),dim(getDimension());
+  ret->alloc(nbGaussPt,nbPtsRefCell*dim);
+  double *retPtr(ret->getPointer());
+  for(int iGaussPt = 0 ; iGaussPt < nbGaussPt ; ++iGaussPt)
+  {
+    std::vector<double> curGaussPt(_gauss_coord.begin()+iGaussPt*dim,_gauss_coord.begin()+(iGaussPt+1)*dim);
+    INTERP_KERNEL::GaussInfo gi(_type,curGaussPt,1,_ref_coord,nbPtsRefCell);
+    gi.initLocalInfo();
+    const double *devOfFuncVal( gi.getDerivativeOfShapeFunctionAt(0) );
+    std::copy(devOfFuncVal,devOfFuncVal+nbPtsRefCell*dim,retPtr);
+    retPtr += nbPtsRefCell*dim;
+  }
+  return ret;
+}
+
 /*!
  * This method sets the comp_th component of ptIdInCell_th point coordinate of reference element of type this->_type.
  * @throw if not 0<=ptIdInCell<nbOfNodePerCell or if not 0<=comp<dim
@@ -324,6 +357,20 @@ bool MEDCouplingGaussLocalization::AreAlmostEqual(const std::vector<double>& v1,
     return false;
   std::vector<double> tmp(sz);
   std::transform(v1.begin(),v1.end(),v2.begin(),tmp.begin(),std::minus<double>());
-  std::transform(tmp.begin(),tmp.end(),tmp.begin(),std::ptr_fun<double,double>(fabs));
+  std::transform(tmp.begin(),tmp.end(),tmp.begin(),[](double c){return fabs(c);});
   return *std::max_element(tmp.begin(),tmp.end())<eps;
 }
+
+MCAuto<DataArrayDouble> MEDCouplingGaussLocalization::GetDefaultReferenceCoordinatesOf(INTERP_KERNEL::NormalizedCellType type)
+{
+  std::vector<double> retCpp(INTERP_KERNEL::GaussInfo::GetDefaultReferenceCoordinatesOf(type));
+  const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
+  auto nbDim(cm.getDimension());
+  std::size_t sz(retCpp.size());
+  MCAuto<DataArrayDouble> ret(DataArrayDouble::New());
+  if( sz%std::size_t(nbDim) != 0 )
+    THROW_IK_EXCEPTION("GetDefaultReferenceCoordinatesOf : unexpected size of defaut array : " << sz << " % " << nbDim << " != 0 !");
+  ret->alloc(sz/size_t(nbDim),nbDim);
+  std::copy(retCpp.begin(),retCpp.end(),ret->getPointer());
+  return ret;
+}