From 5293fda493e627b9815b6c0b0322a4ab828b23c4 Mon Sep 17 00:00:00 2001
From: ageay <ageay>
Date: Fri, 19 Feb 2010 17:27:27 +0000
Subject: [PATCH] Addition of some usefull methods.

---
 src/MEDCoupling/MEDCouplingCMesh.cxx          | 24 +++++++
 src/MEDCoupling/MEDCouplingCMesh.hxx          |  2 +
 src/MEDCoupling/MEDCouplingExtrudedMesh.cxx   | 12 ++++
 src/MEDCoupling/MEDCouplingExtrudedMesh.hxx   |  2 +
 .../MEDCouplingFieldDiscretization.cxx        |  4 +-
 src/MEDCoupling/MEDCouplingFieldDouble.cxx    | 30 +++++++++
 src/MEDCoupling/MEDCouplingFieldDouble.hxx    |  2 +
 src/MEDCoupling/MEDCouplingMesh.hxx           |  2 +
 src/MEDCoupling/MEDCouplingUMesh.cxx          | 64 +++++++++++++++++++
 src/MEDCoupling/MEDCouplingUMesh.hxx          |  2 +
 src/MEDCoupling/MEDCouplingUMeshDesc.cxx      | 14 ++++
 src/MEDCoupling/MEDCouplingUMeshDesc.hxx      |  2 +
 .../Test/MEDCouplingBasicsTest.cxx            | 21 ++++++
 .../Test/MEDCouplingBasicsTest.hxx            |  2 +
 14 files changed, 180 insertions(+), 3 deletions(-)

diff --git a/src/MEDCoupling/MEDCouplingCMesh.cxx b/src/MEDCoupling/MEDCouplingCMesh.cxx
index 34c3f3c8e..523fb31c8 100644
--- a/src/MEDCoupling/MEDCouplingCMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingCMesh.cxx
@@ -18,6 +18,7 @@
 //
 #include "MEDCouplingCMesh.hxx"
 #include "MEDCouplingMemArray.hxx"
+#include "MEDCouplingFieldDouble.hxx"
 
 using namespace ParaMEDMEM;
 
@@ -174,6 +175,29 @@ MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
   return 0;
 }
 
+MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureFieldOnNode(bool isAbs) const
+{
+  //not implemented yet !
+  return 0;
+}
+
+MEDCouplingFieldDouble *MEDCouplingCMesh::buildOrthogonalField() const
+{
+  if(getMeshDimension()!=2)
+    throw INTERP_KERNEL::Exception("Expected a cmesh with meshDim == 2 !");
+  MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
+  DataArrayDouble *array=DataArrayDouble::New();
+  int nbOfCells=getNumberOfCells();
+  array->alloc(nbOfCells,3);
+  double *vals=array->getPointer();
+  for(int i=0;i<nbOfCells;i++)
+    { vals[3*i]=1.; vals[3*i+1]=1.; vals[3*i+2]=1.; }
+  ret->setArray(array);
+  array->decrRef();
+  ret->setMesh(this);
+  return ret;
+}
+
 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
 {
   throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
diff --git a/src/MEDCoupling/MEDCouplingCMesh.hxx b/src/MEDCoupling/MEDCouplingCMesh.hxx
index 7d4f7e03e..4f4096a44 100644
--- a/src/MEDCoupling/MEDCouplingCMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingCMesh.hxx
@@ -46,6 +46,8 @@ namespace ParaMEDMEM
     // tools
     void getBoundingBox(double *bbox) const;
     MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
+    MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
+    MEDCouplingFieldDouble *buildOrthogonalField() const;
     void rotate(const double *center, const double *vector, double angle);
     void translate(const double *vector);
     MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
index 87a61f7c3..f89a0e5d9 100644
--- a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
@@ -120,6 +120,18 @@ MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::getMeasureField(bool) const
   return 0;
 }
 
+MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::getMeasureFieldOnNode(bool) const
+{
+  //not implemented yet
+  return 0;
+}
+
+MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::buildOrthogonalField() const
+{
+  //not implemented yet
+  throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::buildOrthogonalField not implemented yet !");
+}
+
 MEDCouplingExtrudedMesh::~MEDCouplingExtrudedMesh()
 {
   if(_mesh2D)
diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
index 5b6f0ede5..8480bcf2d 100644
--- a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
@@ -47,6 +47,8 @@ namespace ParaMEDMEM
     MEDCouplingUMesh *getMesh1D() const { return _mesh1D; }
     DataArrayInt *getMesh3DIds() const { return _mesh3D_ids; }
     MEDCouplingFieldDouble *getMeasureField(bool) const;
+    MEDCouplingFieldDouble *getMeasureFieldOnNode(bool) const;
+    MEDCouplingFieldDouble *buildOrthogonalField() const;
     static int findCorrespCellByNodalConn(const std::vector<int>& nodalConnec,
                                           const int *revNodalPtr, const int *revNodalIndxPtr) throw(INTERP_KERNEL::Exception);
     void rotate(const double *center, const double *vector, double angle);
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index a2e0493a2..a9f648117 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -179,9 +179,7 @@ void MEDCouplingFieldDiscretizationP1::checkCoherencyBetween(const MEDCouplingMe
 
 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getWeightingField(const MEDCouplingMesh *mesh, bool isAbs) const
 {
-  //not implemented yet.
-  //Dual mesh to build
-  return 0;
+  return mesh->getMeasureFieldOnNode(isAbs);
 }
 
 void MEDCouplingFieldDiscretizationP1::renumberValuesOnNodes(const DataArrayInt *old2New, DataArrayDouble *arr) const
diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
index 5ab8f7998..93beaa328 100644
--- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
@@ -118,6 +118,16 @@ double MEDCouplingFieldDouble::accumulate(int compId) const
   return ret;
 }
 
+void MEDCouplingFieldDouble::accumulate(double *res) const
+{
+  const double *ptr=getArray()->getConstPointer();
+  int nbTuple=getArray()->getNumberOfTuples();
+  int nbComps=getArray()->getNumberOfComponents();
+  std::fill(res,res+nbComps,0.);
+  for(int i=0;i<nbTuple;i++)
+    std::transform(ptr+i*nbComps,ptr+(i+1)*nbComps,res,res,std::plus<double>());
+}
+
 double MEDCouplingFieldDouble::measureAccumulate(int compId, bool isWAbs) const
 {
   if(!_mesh)
@@ -132,6 +142,26 @@ double MEDCouplingFieldDouble::measureAccumulate(int compId, bool isWAbs) const
   return ret;
 }
 
+void MEDCouplingFieldDouble::measureAccumulate(bool isWAbs, double *res) const
+{
+  if(!_mesh)
+    throw INTERP_KERNEL::Exception("No mesh underlying this field to perform measureAccumulate2");
+  MEDCouplingFieldDouble *weight=_type->getWeightingField(_mesh,isWAbs);
+  const double *ptr=weight->getArray()->getConstPointer();
+  int nbOfValues=weight->getArray()->getNbOfElems();
+  int nbComps=getArray()->getNumberOfComponents();
+  const double *vals=getArray()->getConstPointer();
+  std::fill(res,res+nbComps,0.);
+  double *tmp=new double[nbComps];
+  for (int i=0; i<nbOfValues; i++)
+    {
+      std::transform(vals+i*nbComps,vals+(i+1)*nbComps,tmp,std::bind2nd(std::multiplies<double>(),ptr[i]));
+      std::transform(tmp,tmp+nbComps,res,res,std::plus<double>());
+    }
+  weight->decrRef();
+  delete [] tmp;
+}
+
 void MEDCouplingFieldDouble::getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception)
 {
   _time_discr->checkNoTimePresence();
diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx
index 16496f0ec..bcece9122 100644
--- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx
@@ -49,7 +49,9 @@ namespace ParaMEDMEM
     void setArray(DataArrayDouble *array);
     DataArrayDouble *getArray() const { return _time_discr->getArray(); }
     double accumulate(int compId) const;
+    void accumulate(double *res) const;
     double measureAccumulate(int compId, bool isWAbs) const;
+    void measureAccumulate(bool isWAbs, double *res) const;
     void getValueOn(const double *spaceLoc, double *res) const throw(INTERP_KERNEL::Exception);
     void getValueOn(const double *spaceLoc, double time, double *res) const throw(INTERP_KERNEL::Exception);
     //! \b temporary
diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx
index f49496d6d..e28781827 100644
--- a/src/MEDCoupling/MEDCouplingMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingMesh.hxx
@@ -55,7 +55,9 @@ namespace ParaMEDMEM
     // tools
     virtual void getBoundingBox(double *bbox) const = 0;
     virtual MEDCouplingFieldDouble *getMeasureField(bool isAbs) const = 0;
+    virtual MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const = 0;
     virtual MEDCouplingFieldDouble *fillFromAnalytic(TypeOfField t, int nbOfComp, FunctionToEvaluate func) const;
+    virtual MEDCouplingFieldDouble *buildOrthogonalField() const = 0;
     virtual void rotate(const double *center, const double *vector, double angle) = 0;
     virtual void translate(const double *vector) = 0;
     virtual MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const = 0;
diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx
index 984704e1a..5bb6c99e2 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.cxx
@@ -20,6 +20,7 @@
 #include "MEDCouplingFieldDouble.hxx"
 #include "CellModel.hxx"
 #include "VolSurfUser.txx"
+#include "InterpolationUtils.hxx"
 
 #include <sstream>
 #include <limits>
@@ -840,6 +841,69 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField(bool isAbs) const
   return field;
 }
 
+MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureFieldOnNode(bool isAbs) const
+{
+  MEDCouplingFieldDouble *tmp=getMeasureField(abs);
+  std::string name="MeasureOnNodeOfMesh_";
+  name+=getName();
+  int nbNodes=getNumberOfNodes();
+  MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_NODES);
+  double cst=1./((double)getMeshDimension()+1.);
+  DataArrayDouble* array=DataArrayDouble::New();
+  array->alloc(nbNodes,1);
+  double *valsToFill=array->getPointer();
+  std::fill(valsToFill,valsToFill+nbNodes,0.);
+  const double *values=tmp->getArray()->getConstPointer();
+  DataArrayInt *da=DataArrayInt::New();
+  DataArrayInt *daInd=DataArrayInt::New();
+  getReverseNodalConnectivity(da,daInd);
+  const int *daPtr=da->getConstPointer();
+  const int *daIPtr=daInd->getConstPointer();
+  for(int i=0;i<nbNodes;i++)
+    for(const int *cell=daPtr+daIPtr[i];cell!=daPtr+daIPtr[i+1];cell++)
+      valsToFill[i]+=cst*values[*cell];
+  ret->setMesh(this);
+  da->decrRef();
+  daInd->decrRef();
+  ret->setArray(array);
+  array->decrRef();
+  tmp->decrRef();
+  return ret;
+}
+
+MEDCouplingFieldDouble *MEDCouplingUMesh::buildOrthogonalField() const
+{
+  if(getMeshDimension()!=2)
+    throw INTERP_KERNEL::Exception("Expected a umesh with meshDim == 2 !");
+  MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
+  DataArrayDouble *array=DataArrayDouble::New();
+  int nbOfCells=getNumberOfCells();
+  array->alloc(nbOfCells,3);
+  double *vals=array->getPointer();
+  const int *connI=_nodal_connec_index->getConstPointer();
+  const int *conn=_nodal_connec->getConstPointer();
+  const double *coords=_coords->getConstPointer();
+  if(getSpaceDimension()==3)
+    {
+      for(int i=0;i<nbOfCells;i++,vals+=3)
+        {
+          int offset=connI[i];
+          INTERP_KERNEL::crossprod<3>(coords+3*conn[offset+1],coords+3*conn[offset+2],coords+3*conn[offset+3],vals);
+          double n=INTERP_KERNEL::norm<3>(vals);
+          std::transform(vals,vals+3,vals,std::bind2nd(std::multiplies<double>(),1./n));
+        }
+    }
+  else
+    {
+      for(int i=0;i<nbOfCells;i++)
+        { vals[3*i]=1.; vals[3*i+1]=1.; vals[3*i+2]=1.; }
+    }
+  ret->setArray(array);
+  array->decrRef();
+  ret->setMesh(this);
+  return ret;
+}
+
 /*!
  * This method is only available for a mesh with meshDim==2 and spaceDim==2||spaceDim==3.
  * This method returns a vector 'cells' where all detected butterfly cells have been added to cells.
diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx
index e793c931c..9c036342f 100644
--- a/src/MEDCoupling/MEDCouplingUMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingUMesh.hxx
@@ -70,6 +70,8 @@ namespace ParaMEDMEM
     void renumberNodes(const int *newNodeNumbers, int newNbOfNodes);
     void giveElemsInBoundingBox(const double *bbox, double eps, std::vector<int>& elems);
     MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
+    MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
+    MEDCouplingFieldDouble *buildOrthogonalField() const;
     void checkButterflyCells(std::vector<int>& cells) const;
     bool checkConsecutiveCellTypes() const;
     MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
index 1dba15e9a..3a34d1f84 100644
--- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
+++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
@@ -262,6 +262,20 @@ MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureField(bool isAbs) const
   return 0;
 }
 
+MEDCouplingFieldDouble *MEDCouplingUMeshDesc::getMeasureFieldOnNode(bool isAbs) const
+{
+  //not implemented yet.
+  return 0;
+}
+
+MEDCouplingFieldDouble *MEDCouplingUMeshDesc::buildOrthogonalField() const
+{
+  if(getMeshDimension()!=2)
+    throw INTERP_KERNEL::Exception("Expected a cmesh with meshDim == 2 !");
+  //not implemented yet !
+  return 0;
+}
+
 DataArrayInt *MEDCouplingUMeshDesc::zipCoordsTraducer()
 {
   //not implemented yet.
diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
index b6ef8aa0c..ebd34115a 100644
--- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
+++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
@@ -56,6 +56,8 @@ namespace ParaMEDMEM
     MEDCouplingPointSet *buildBoundaryMesh(bool keepCoords) const;
     void renumberNodes(const int *newNodeNumbers, int newNbOfNodes);
     MEDCouplingFieldDouble *getMeasureField(bool isAbs) const;
+    MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
+    MEDCouplingFieldDouble *buildOrthogonalField() const;
     DataArrayInt *zipCoordsTraducer();
     MEDCouplingMesh *mergeMyselfWith(const MEDCouplingMesh *other) const;
     DataArrayDouble *getBarycenterAndOwner() const;
diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
index 5b5263c12..bfa87a1a3 100644
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.cxx
@@ -1140,6 +1140,13 @@ void MEDCouplingBasicsTest::testFillFromAnalytic()
   std::transform(values3,values3+18,values3,std::ptr_fun<double,double>(fabs));
   max=*std::max_element(values3,values3+18);
   CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,max,1.e-12);
+  double values4[2];
+  f1->accumulate(values4);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(3.6,values4[0],1.e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(7.2,values4[1],1.e-12);
+  f1->measureAccumulate(true,values4);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(0.5,values4[0],1.e-12);
+  CPPUNIT_ASSERT_DOUBLES_EQUAL(1.,values4[1],1.e-12);
   f1->decrRef();
   //
   CPPUNIT_ASSERT_THROW(f1=m->fillFromAnalytic(ON_NODES,1,func3),INTERP_KERNEL::Exception);
@@ -1317,6 +1324,20 @@ void MEDCouplingBasicsTest::testCheckConsecutiveCellTypes()
   sourceMesh->decrRef();
 }
 
+void MEDCouplingBasicsTest::testBuildOrthogonalField()
+{
+  MEDCouplingUMesh *targetMesh=build3DSurfTargetMesh_1();
+  MEDCouplingFieldDouble *field=targetMesh->buildOrthogonalField();
+  double expected[3]={0.70710678118654746,0.,-0.70710678118654746};
+  CPPUNIT_ASSERT_EQUAL(5,field->getNumberOfTuples());
+  CPPUNIT_ASSERT_EQUAL(3,field->getNumberOfComponents());
+  const double *vals=field->getArray()->getConstPointer();
+  for(int i=0;i<15;i++)
+    CPPUNIT_ASSERT_DOUBLES_EQUAL(expected[i%3],vals[i],1e-12);
+  field->decrRef();
+  targetMesh->decrRef();
+}
+
 void MEDCouplingBasicsTest::test2DInterpP0P0_1()
 {
   MEDCouplingUMesh *sourceMesh=build2DSourceMesh_1();
diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
index ef9558d4b..9d3628313 100644
--- a/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
+++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest.hxx
@@ -59,6 +59,7 @@ namespace ParaMEDMEM
     CPPUNIT_TEST( testOperationsOnFields );
     CPPUNIT_TEST( testMergeNodesOnField );
     CPPUNIT_TEST( testCheckConsecutiveCellTypes );
+    CPPUNIT_TEST( testBuildOrthogonalField );
     CPPUNIT_TEST( test2DInterpP0P0_1 );
     CPPUNIT_TEST( test2DInterpP0P0PL_1 );
     CPPUNIT_TEST( test2DInterpP0P0PL_2 );
@@ -137,6 +138,7 @@ namespace ParaMEDMEM
     void testOperationsOnFields();
     void testMergeNodesOnField();
     void testCheckConsecutiveCellTypes();
+    void testBuildOrthogonalField();
     void test2DInterpP0P0_1();
     void test2DInterpP0P0PL_1();
     void test2DInterpP0P0PL_2();
-- 
2.39.2