#include "MEDCouplingExtrudedMesh.hxx"
#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingMemArray.hxx"
+#include "MEDCouplingFieldDouble.hxx"
#include "CellModel.hxx"
+#include "InterpolationUtils.hxx"
+
#include <limits>
#include <algorithm>
#include <functional>
void MEDCouplingExtrudedMesh::project1DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
MEDCouplingUMesh *&m1r, MEDCouplingUMesh *&m2r, double *v) throw(INTERP_KERNEL::Exception)
{
+ if(m1->getSpaceDimension()!=3 || m1->getSpaceDimension()!=3)
+ throw INTERP_KERNEL::Exception("Input meshes are expected to have a spaceDim==3 for projec1D !");
m1r=m1->clone(true);
m2r=m2->clone(true);
m1r->changeSpaceDimension(1);
m2r->changeSpaceDimension(1);
+ std::vector<int> c;
+ std::vector<double> ref,ref2;
+ m1->getNodeIdsOfCell(0,c);
+ m1->getCoordinatesOfNode(c[0],ref);
+ m1->getCoordinatesOfNode(c[1],ref2);
+ std::transform(ref2.begin(),ref2.end(),ref.begin(),v,std::minus<double>());
+ double n=INTERP_KERNEL::norm<3>(v);
+ std::transform(v,v+3,v,std::bind2nd(std::multiplies<double>(),1/n));
+ m1->project1D(&ref[0],v,eps,m1r->getCoords()->getPointer());
+ m2->project1D(&ref[0],v,eps,m2r->getCoords()->getPointer());
+
}
void MEDCouplingExtrudedMesh::rotate(const double *center, const double *vector, double angle)
#include "MEDCouplingExtrudedMesh.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
+#include "Interpolation1D.txx"
#include "Interpolation2DCurve.txx"
#include "Interpolation2D.txx"
#include "Interpolation3D.txx"
INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh->getMesh2D());
MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh->getMesh2D());
- INTERP_KERNEL::Interpolation2D interpolation(*this);
+ INTERP_KERNEL::Interpolation2D interpolation2D(*this);
std::vector<std::map<int,double> > matrix2D;
- int nbCols2D=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method);
+ int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method);
MEDCouplingUMesh *s1D,*t1D;
double v[3];
MEDCouplingExtrudedMesh::project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
- MEDCouplingNormalizedUnstructuredMesh<2,1> s1DWrapper(s1D);
- MEDCouplingNormalizedUnstructuredMesh<2,1> t1DWrapper(t1D);
+ MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
+ MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
std::vector<std::map<int,double> > matrix1D;
- int nbCols1D=interpolation.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method);
- INTERP_KERNEL::Interpolation2DCurve myInterpolator;
- //
- _matrix.resize(matrix2D.size()*matrix1D.size());
+ INTERP_KERNEL::Interpolation1D interpolation1D(*this);
+ int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method);
+ s1D->decrRef();
+ t1D->decrRef();
+ buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
+ target_mesh->getMesh3DIds()->getConstPointer());
//
_deno_multiply.clear();
_deno_multiply.resize(_matrix.size());
deno[idx][(*iter2).first]=values[(*iter2).first];
}
}
+
+void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
+ const std::vector< std::map<int,double> >& m2D,
+ const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
+ const int *corrCellIdTrg)
+{
+ int nbOf2DCellsTrg=m2D.size();
+ int nbOf1DCellsTrg=m1D.size();
+ int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
+ _matrix.resize(nbOf3DCellsTrg);
+ int id2R=0;
+ for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
+ {
+ for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
+ {
+ int id1R=0;
+ for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
+ {
+ for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
+ {
+ _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
+ }
+ }
+ }
+ }
+}
void computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField);
void computeProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer);
void computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer);
+ void buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
+ const std::vector< std::map<int,double> >& m2D,
+ const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
+ const int *corrCellIdTrg);
static void reverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn,
std::vector<std::map<int,double> >& matOut);
static void computeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
#include "BBTree.txx"
#include <sstream>
+#include <numeric>
#include <limits>
#include <list>
getNodeIdsOfCell(i,conn);
pt=std::transform(coo+conn[1]*spaceDim,coo+(conn[1]+1)*spaceDim,coo+conn[0]*spaceDim,pt,std::minus<double>());
}
- array->alloc(nbOfCells,spaceDim);
+ ret->setArray(array);
array->decrRef();
ret->setMesh(this);
return ret;
}
+/*!
+ * This method is only callable on mesh with meshdim == 1 containing only SEG2 and spaceDim==3.
+ * This method projects this on the 3D line defined by (pt,v). This methods first checks that all SEG2 are along v vector.
+ * @param pt reference point of the line
+ * @param v normalized director vector of the line
+ * @param eps max precision before throwing an exception
+ * @param res output of size this->getNumberOfCells
+ */
+void MEDCouplingUMesh::project1D(const double *pt, const double *v, double eps, double *res) const
+{
+ if(getMeshDimension()!=1)
+ throw INTERP_KERNEL::Exception("Expected a umesh with meshDim == 1 for project1D !");
+ if(_types.size()!=1 || *(_types.begin())!=INTERP_KERNEL::NORM_SEG2)
+ throw INTERP_KERNEL::Exception("Expected a umesh with only NORM_SEG2 type of elements for project1D !");
+ if(getSpaceDimension()!=3)
+ throw INTERP_KERNEL::Exception("Expected a umesh with spaceDim==3 for project1D !");
+ MEDCouplingFieldDouble *f=buildLinearField();
+ const double *fPtr=f->getArray()->getConstPointer();
+ double tmp[3];
+ for(int i=0;i<getNumberOfCells();i++)
+ {
+ const double *tmp1=fPtr+3*i;
+ tmp[0]=tmp1[1]*v[2]-tmp1[2]*v[1];
+ tmp[1]=tmp1[2]*v[0]-tmp1[0]*v[2];
+ tmp[2]=tmp1[0]*v[1]-tmp1[1]*v[0];
+ double n1=INTERP_KERNEL::norm<3>(tmp);
+ n1/=INTERP_KERNEL::norm<3>(tmp1);
+ if(n1>eps)
+ {
+ f->decrRef();
+ throw INTERP_KERNEL::Exception("UMesh::Projection 1D failed !");
+ }
+ }
+ const double *coo=getCoords()->getConstPointer();
+ for(int i=0;i<getNumberOfNodes();i++)
+ {
+ std::transform(coo+i*3,coo+i*3+3,pt,tmp,std::minus<double>());
+ std::transform(tmp,tmp+3,v,tmp,std::multiplies<double>());
+ res[i]=std::accumulate(tmp,tmp+3,0.);
+ }
+ f->decrRef();
+}
+
/*!
* Returns a cell if any that contains the point located on 'pos' with precison eps.
* If 'pos' is outside 'this' -1 is returned. If several cells contain this point the cell with the smallest id is returned.
MEDCOUPLING_EXPORT MEDCouplingFieldDouble *getMeasureFieldOnNode(bool isAbs) const;
MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildOrthogonalField() const;
MEDCOUPLING_EXPORT MEDCouplingFieldDouble *buildLinearField() const;
+ MEDCOUPLING_EXPORT void project1D(const double *pt, const double *v, double eps, double *res) const;
MEDCOUPLING_EXPORT int getCellContainingPoint(const double *pos, double eps) const;
MEDCOUPLING_EXPORT void getCellsContainingPoint(const double *pos, double eps, std::vector<int>& elts) const;
MEDCOUPLING_EXPORT void getCellsContainingPoints(const double *pos, int nbOfPoints, double eps, std::vector<int>& elts, std::vector<int>& eltsIndex) const;
