]> SALOME platform Git repositories - tools/medcoupling.git/blob - src/MEDCoupling/MEDCouplingUMesh_internal.hxx
Salome HOME
Refactoring QuadraticPlanar precision tuning. Fixes intersector bugs.
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh_internal.hxx
1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // Author : Anthony Geay (EdF)
20
21 using namespace MEDCoupling;
22
23 class MinusOneSonsGenerator
24 {
25 public:
26   MinusOneSonsGenerator(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
27   unsigned getNumberOfSons2(const int *conn, int lgth) const { return _cm.getNumberOfSons2(conn,lgth); }
28   unsigned fillSonCellNodalConnectivity2(int sonId, const int *nodalConn, int lgth, int *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonCellNodalConnectivity2(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
29   static const int DELTA=1;
30 private:
31   const INTERP_KERNEL::CellModel& _cm;
32 };
33
34 class MinusOneSonsGeneratorBiQuadratic
35 {
36 public:
37   MinusOneSonsGeneratorBiQuadratic(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
38   unsigned getNumberOfSons2(const int *conn, int lgth) const { return _cm.getNumberOfSons2(conn,lgth); }
39   unsigned fillSonCellNodalConnectivity2(int sonId, const int *nodalConn, int lgth, int *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonCellNodalConnectivity4(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
40   static const int DELTA=1;
41 private:
42   const INTERP_KERNEL::CellModel& _cm;
43 };
44
45 class MinusTwoSonsGenerator
46 {
47 public:
48   MinusTwoSonsGenerator(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
49   unsigned getNumberOfSons2(const int *conn, int lgth) const { return _cm.getNumberOfEdgesIn3D(conn,lgth); }
50   unsigned fillSonCellNodalConnectivity2(int sonId, const int *nodalConn, int lgth, int *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonEdgesNodalConnectivity3D(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
51   static const int DELTA=2;
52 private:
53   const INTERP_KERNEL::CellModel& _cm;
54 };
55
56 class MicroEdgesGenerator2D
57 {
58 public:
59   MicroEdgesGenerator2D(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
60   unsigned getNumberOfSons2(const int *conn, int lgth) const { return _cm.getNumberOfMicroEdges(); }
61   unsigned fillSonCellNodalConnectivity2(int sonId, const int *nodalConn, int lgth, int *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillMicroEdgeNodalConnectivity(sonId,nodalConn,sonNodalConn,typeOfSon); }
62   static const int DELTA=1;
63 private:
64   const INTERP_KERNEL::CellModel& _cm;
65 };
66
67 class MicroEdgesGenerator3D
68 {
69 public:
70   MicroEdgesGenerator3D(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
71   unsigned getNumberOfSons2(const int *conn, int lgth) const { return _cm.getNumberOfMicroEdges(); }
72   unsigned fillSonCellNodalConnectivity2(int sonId, const int *nodalConn, int lgth, int *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillMicroEdgeNodalConnectivity(sonId,nodalConn,sonNodalConn,typeOfSon); }
73   static const int DELTA=2;
74 private:
75   const INTERP_KERNEL::CellModel& _cm;
76 };
77
78 int MEDCouplingFastNbrer(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2);
79 int MEDCouplingOrientationSensitiveNbrer(int id, unsigned nb, const INTERP_KERNEL::CellModel& cm, bool compute, const int *conn1, const int *conn2);
80
81 namespace MEDCoupling
82 {
83   template<const int SPACEDIMM>
84   class DummyClsMCUG
85   {
86   public:
87     static const int MY_SPACEDIM=SPACEDIMM;
88     static const int MY_MESHDIM=8;
89     typedef int MyConnType;
90     static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
91     // begin
92     // useless, but for windows compilation ...
93     const double* getCoordinatesPtr() const { return 0; }
94     const int* getConnectivityPtr() const { return 0; }
95     const int* getConnectivityIndexPtr() const { return 0; }
96     INTERP_KERNEL::NormalizedCellType getTypeOfElement(int) const { return (INTERP_KERNEL::NormalizedCellType)0; }
97     // end
98   };
99 }
100
101 template<int SPACEDIM>
102 void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const double *pos, int nbOfPoints,
103                                                    double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const
104 {
105   // Override precision for this method only:
106   INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
107   INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(eps);
108
109   elts=DataArrayInt::New(); eltsIndex=DataArrayInt::New(); eltsIndex->alloc(nbOfPoints+1,1); eltsIndex->setIJ(0,0,0); elts->alloc(0,1);
110   int *eltsIndexPtr(eltsIndex->getPointer());
111   MCAuto<DataArrayDouble> bboxArr(getBoundingBoxForBBTree(eps));
112   const double *bbox(bboxArr->begin());
113   int nbOfCells=getNumberOfCells();
114   const int *conn=_nodal_connec->getConstPointer();
115   const int *connI=_nodal_connec_index->getConstPointer();
116   double bb[2*SPACEDIM];
117   BBTree<SPACEDIM,int> myTree(&bbox[0],0,0,nbOfCells,-eps);
118   for(int i=0;i<nbOfPoints;i++)
119     {
120       eltsIndexPtr[i+1]=eltsIndexPtr[i];
121       for(int j=0;j<SPACEDIM;j++)
122         {
123           bb[2*j]=pos[SPACEDIM*i+j];
124           bb[2*j+1]=pos[SPACEDIM*i+j];
125         }
126       std::vector<int> candidates;
127       myTree.getIntersectingElems(bb,candidates);
128       for(std::vector<int>::const_iterator iter=candidates.begin();iter!=candidates.end();iter++)
129         {
130           int sz(connI[(*iter)+1]-connI[*iter]-1);
131           INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]]);
132           bool status(false);
133           if(ct!=INTERP_KERNEL::NORM_POLYGON && ct!=INTERP_KERNEL::NORM_QPOLYG)
134             status=INTERP_KERNEL::PointLocatorAlgos<DummyClsMCUG<SPACEDIM> >::isElementContainsPoint(pos+i*SPACEDIM,ct,coords,conn+connI[*iter]+1,sz,eps);
135           else
136             {
137               if(SPACEDIM!=2)
138                 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPointsAlg : not implemented yet for POLYGON and QPOLYGON in spaceDim 3 !");
139               std::vector<INTERP_KERNEL::Node *> nodes(sz);
140               INTERP_KERNEL::QuadraticPolygon *pol(0);
141               for(int j=0;j<sz;j++)
142                 {
143                   int nodeId(conn[connI[*iter]+1+j]);
144                   nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*SPACEDIM],coords[nodeId*SPACEDIM+1]);
145                 }
146               if(!INTERP_KERNEL::CellModel::GetCellModel(ct).isQuadratic())
147                 pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
148               else
149                 pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
150               INTERP_KERNEL::Node *n(new INTERP_KERNEL::Node(pos[i*SPACEDIM],pos[i*SPACEDIM+1]));
151               double a(0.),b(0.),c(0.);
152               a=pol->normalizeMe(b,c); n->applySimilarity(b,c,a);
153               status=pol->isInOrOut2(n);
154               delete pol; n->decrRef();
155             }
156           if(status)
157             {
158               eltsIndexPtr[i+1]++;
159               elts->pushBackSilent(*iter);
160             }
161         }
162     }
163 }
164
165 /*!
166  * \b WARNING this method do the assumption that connectivity lies on the coordinates set.
167  * For speed reasons no check of this will be done.
168  */
169 template<class SonsGenerator>
170 MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivityGen(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx, DimM1DescNbrer nbrer) const
171 {
172   if(!desc || !descIndx || !revDesc || !revDescIndx)
173     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildDescendingConnectivityGen : present of a null pointer in input !");
174   checkConnectivityFullyDefined();
175   int nbOfCells=getNumberOfCells();
176   int nbOfNodes=getNumberOfNodes();
177   MCAuto<DataArrayInt> revNodalIndx=DataArrayInt::New(); revNodalIndx->alloc(nbOfNodes+1,1); revNodalIndx->fillWithZero();
178   int *revNodalIndxPtr=revNodalIndx->getPointer();
179   const int *conn=_nodal_connec->getConstPointer();
180   const int *connIndex=_nodal_connec_index->getConstPointer();
181   std::string name="Mesh constituent of "; name+=getName();
182   MCAuto<MEDCouplingUMesh> ret=MEDCouplingUMesh::New(name,getMeshDimension()-SonsGenerator::DELTA);
183   ret->setCoords(getCoords());
184   ret->allocateCells(2*nbOfCells);
185   descIndx->alloc(nbOfCells+1,1);
186   MCAuto<DataArrayInt> revDesc2(DataArrayInt::New()); revDesc2->reserve(2*nbOfCells);
187   int *descIndxPtr=descIndx->getPointer(); *descIndxPtr++=0;
188   for(int eltId=0;eltId<nbOfCells;eltId++,descIndxPtr++)
189     {
190       int pos=connIndex[eltId];
191       int posP1=connIndex[eltId+1];
192       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[pos]);
193       SonsGenerator sg(cm);
194       unsigned nbOfSons=sg.getNumberOfSons2(conn+pos+1,posP1-pos-1);
195       INTERP_KERNEL::AutoPtr<int> tmp=new int[posP1-pos];
196       for(unsigned i=0;i<nbOfSons;i++)
197         {
198           INTERP_KERNEL::NormalizedCellType cmsId;
199           unsigned nbOfNodesSon=sg.fillSonCellNodalConnectivity2(i,conn+pos+1,posP1-pos-1,tmp,cmsId);
200           for(unsigned k=0;k<nbOfNodesSon;k++)
201             if(tmp[k]>=0)
202               revNodalIndxPtr[tmp[k]+1]++;
203           ret->insertNextCell(cmsId,nbOfNodesSon,tmp);
204           revDesc2->pushBackSilent(eltId);
205         }
206       descIndxPtr[0]=descIndxPtr[-1]+(int)nbOfSons;
207     }
208   int nbOfCellsM1=ret->getNumberOfCells();
209   std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<int>());
210   MCAuto<DataArrayInt> revNodal=DataArrayInt::New(); revNodal->alloc(revNodalIndx->back(),1);
211   std::fill(revNodal->getPointer(),revNodal->getPointer()+revNodalIndx->back(),-1);
212   int *revNodalPtr=revNodal->getPointer();
213   const int *connM1=ret->getNodalConnectivity()->getConstPointer();
214   const int *connIndexM1=ret->getNodalConnectivityIndex()->getConstPointer();
215   for(int eltId=0;eltId<nbOfCellsM1;eltId++)
216     {
217       const int *strtNdlConnOfCurCell=connM1+connIndexM1[eltId]+1;
218       const int *endNdlConnOfCurCell=connM1+connIndexM1[eltId+1];
219       for(const int *iter=strtNdlConnOfCurCell;iter!=endNdlConnOfCurCell;iter++)
220         if(*iter>=0)//for polyhedrons
221           *std::find_if(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1],std::bind2nd(std::equal_to<int>(),-1))=eltId;
222     }
223   //
224   DataArrayInt *commonCells=0,*commonCellsI=0;
225   FindCommonCellsAlg(3,0,ret->getNodalConnectivity(),ret->getNodalConnectivityIndex(),revNodal,revNodalIndx,commonCells,commonCellsI);
226   MCAuto<DataArrayInt> commonCellsTmp(commonCells),commonCellsITmp(commonCellsI);
227   const int *commonCellsPtr(commonCells->getConstPointer()),*commonCellsIPtr(commonCellsI->getConstPointer());
228   int newNbOfCellsM1=-1;
229   MCAuto<DataArrayInt> o2nM1=DataArrayInt::ConvertIndexArrayToO2N(nbOfCellsM1,commonCells->begin(),
230                                                                                                             commonCellsI->begin(),commonCellsI->end(),newNbOfCellsM1);
231   std::vector<bool> isImpacted(nbOfCellsM1,false);
232   for(const int *work=commonCellsI->begin();work!=commonCellsI->end()-1;work++)
233     for(int work2=work[0];work2!=work[1];work2++)
234       isImpacted[commonCellsPtr[work2]]=true;
235   const int *o2nM1Ptr=o2nM1->getConstPointer();
236   MCAuto<DataArrayInt> n2oM1=o2nM1->invertArrayO2N2N2OBis(newNbOfCellsM1);
237   const int *n2oM1Ptr=n2oM1->getConstPointer();
238   MCAuto<MEDCouplingUMesh> ret2=static_cast<MEDCouplingUMesh *>(ret->buildPartOfMySelf(n2oM1->begin(),n2oM1->end(),true));
239   ret2->copyTinyInfoFrom(this);
240   desc->alloc(descIndx->back(),1);
241   int *descPtr=desc->getPointer();
242   const INTERP_KERNEL::CellModel& cmsDft=INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_POINT1);
243   for(int i=0;i<nbOfCellsM1;i++,descPtr++)
244     {
245       if(!isImpacted[i])
246         *descPtr=nbrer(o2nM1Ptr[i],0,cmsDft,false,0,0);
247       else
248         {
249           if(i!=n2oM1Ptr[o2nM1Ptr[i]])
250             {
251               const INTERP_KERNEL::CellModel& cms=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)connM1[connIndexM1[i]]);
252               *descPtr=nbrer(o2nM1Ptr[i],connIndexM1[i+1]-connIndexM1[i]-1,cms,true,connM1+connIndexM1[n2oM1Ptr[o2nM1Ptr[i]]]+1,connM1+connIndexM1[i]+1);
253             }
254           else
255             *descPtr=nbrer(o2nM1Ptr[i],0,cmsDft,false,0,0);
256         }
257     }
258   revDesc->reserve(newNbOfCellsM1);
259   revDescIndx->alloc(newNbOfCellsM1+1,1);
260   int *revDescIndxPtr=revDescIndx->getPointer(); *revDescIndxPtr++=0;
261   const int *revDesc2Ptr=revDesc2->getConstPointer();
262   for(int i=0;i<newNbOfCellsM1;i++,revDescIndxPtr++)
263     {
264       int oldCellIdM1=n2oM1Ptr[i];
265       if(!isImpacted[oldCellIdM1])
266         {
267           revDesc->pushBackSilent(revDesc2Ptr[oldCellIdM1]);
268           revDescIndxPtr[0]=revDescIndxPtr[-1]+1;
269         }
270       else
271         {
272           for(int j=commonCellsIPtr[0];j<commonCellsIPtr[1];j++)
273             revDesc->pushBackSilent(revDesc2Ptr[commonCellsPtr[j]]);
274           revDescIndxPtr[0]=revDescIndxPtr[-1]+commonCellsIPtr[1]-commonCellsIPtr[0];
275           commonCellsIPtr++;
276         }
277     }
278   //
279   return ret2.retn();
280 }
281
282