1 // Copyright (C) 2007-2020 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (EdF)
21 using namespace MEDCoupling;
23 class MinusOneSonsGenerator
26 MinusOneSonsGenerator(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
27 unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfSons2(conn,lgth); }
28 unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonCellNodalConnectivity2(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
29 static const int DELTA=1;
31 const INTERP_KERNEL::CellModel& _cm;
34 class MinusOneSonsGeneratorBiQuadratic
37 MinusOneSonsGeneratorBiQuadratic(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
38 unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfSons2(conn,lgth); }
39 unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonCellNodalConnectivity4(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
40 static const int DELTA=1;
42 const INTERP_KERNEL::CellModel& _cm;
45 class MinusTwoSonsGenerator
48 MinusTwoSonsGenerator(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
49 unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfEdgesIn3D(conn,lgth); }
50 unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillSonEdgesNodalConnectivity3D(sonId,nodalConn,lgth,sonNodalConn,typeOfSon); }
51 static const int DELTA=2;
53 const INTERP_KERNEL::CellModel& _cm;
56 class MicroEdgesGenerator2D
59 MicroEdgesGenerator2D(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
60 unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfMicroEdges(); }
61 unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillMicroEdgeNodalConnectivity(sonId,nodalConn,sonNodalConn,typeOfSon); }
62 static const int DELTA=1;
64 const INTERP_KERNEL::CellModel& _cm;
67 class MicroEdgesGenerator3D
70 MicroEdgesGenerator3D(const INTERP_KERNEL::CellModel& cm):_cm(cm) { }
71 unsigned getNumberOfSons2(const mcIdType *conn, mcIdType lgth) const { return _cm.getNumberOfMicroEdges(); }
72 unsigned fillSonCellNodalConnectivity2(int sonId, const mcIdType *nodalConn, mcIdType lgth, mcIdType *sonNodalConn, INTERP_KERNEL::NormalizedCellType& typeOfSon) const { return _cm.fillMicroEdgeNodalConnectivity(sonId,nodalConn,sonNodalConn,typeOfSon); }
73 static const int DELTA=2;
75 const INTERP_KERNEL::CellModel& _cm;
78 mcIdType MEDCouplingFastNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2);
79 mcIdType MEDCouplingOrientationSensitiveNbrer(mcIdType id, mcIdType nb, const INTERP_KERNEL::CellModel& cm, bool compute, const mcIdType *conn1, const mcIdType *conn2);
83 template<const int SPACEDIMM>
87 static const int MY_SPACEDIM=SPACEDIMM;
88 static const int MY_MESHDIM=8;
89 typedef mcIdType MyConnType;
90 static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
92 // useless, but for windows compilation ...
93 const double* getCoordinatesPtr() const { return 0; }
94 const MyConnType* getConnectivityPtr() const { return 0; }
95 const MyConnType* getConnectivityIndexPtr() const { return 0; }
96 INTERP_KERNEL::NormalizedCellType getTypeOfElement(MyConnType) const { return (INTERP_KERNEL::NormalizedCellType)0; }
101 template<int SPACEDIM>
102 void MEDCouplingUMesh::getCellsContainingPointsAlg(const double *coords, const double *pos, mcIdType nbOfPoints,
103 double eps, MCAuto<DataArrayIdType>& elts, MCAuto<DataArrayIdType>& eltsIndex, std::function<bool(INTERP_KERNEL::NormalizedCellType,int)> sensibilityTo2DQuadraticLinearCellsFunc) const
105 // Override precision for this method only:
106 INTERP_KERNEL::QuadraticPlanarPrecision prec(eps);
108 elts=DataArrayIdType::New(); eltsIndex=DataArrayIdType::New(); eltsIndex->alloc(nbOfPoints+1,1); eltsIndex->setIJ(0,0,0); elts->alloc(0,1);
109 mcIdType *eltsIndexPtr(eltsIndex->getPointer());
110 MCAuto<DataArrayDouble> bboxArr(getBoundingBoxForBBTree(eps));
111 const double *bbox(bboxArr->begin());
112 mcIdType nbOfCells=getNumberOfCells();
113 const mcIdType *conn=_nodal_connec->getConstPointer();
114 const mcIdType *connI=_nodal_connec_index->getConstPointer();
115 double bb[2*SPACEDIM];
116 BBTree<SPACEDIM,mcIdType> myTree(&bbox[0],0,0,nbOfCells,-eps);
117 for(mcIdType i=0;i<nbOfPoints;i++)
119 eltsIndexPtr[i+1]=eltsIndexPtr[i];
120 for(int j=0;j<SPACEDIM;j++)
122 bb[2*j]=pos[SPACEDIM*i+j];
123 bb[2*j+1]=pos[SPACEDIM*i+j];
125 std::vector<mcIdType> candidates;
126 myTree.getIntersectingElems(bb,candidates);
127 for(std::vector<mcIdType>::const_iterator iter=candidates.begin();iter!=candidates.end();iter++)
129 mcIdType sz(connI[(*iter)+1]-connI[*iter]-1);
130 INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)conn[connI[*iter]]);
132 // [ABN] : point locator algorithms are only properly working for linear cells.
133 if(ct!=INTERP_KERNEL::NORM_POLYGON && !sensibilityTo2DQuadraticLinearCellsFunc(ct,_mesh_dim))
134 status=INTERP_KERNEL::PointLocatorAlgos<DummyClsMCUG<SPACEDIM> >::isElementContainsPoint(pos+i*SPACEDIM,ct,coords,conn+connI[*iter]+1,sz,eps);
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(mcIdType j=0;j<sz;j++)
143 mcIdType nodeId(conn[connI[*iter]+1+j]);
144 nodes[j]=new INTERP_KERNEL::Node(coords[nodeId*SPACEDIM],coords[nodeId*SPACEDIM+1]);
146 if(!INTERP_KERNEL::CellModel::GetCellModel(ct).isQuadratic())
147 pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes);
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();
159 elts->pushBackSilent(*iter);
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.
169 template<class SonsGenerator>
170 MEDCouplingUMesh *MEDCouplingUMesh::buildDescendingConnectivityGen(DataArrayIdType *desc, DataArrayIdType *descIndx, DataArrayIdType *revDesc, DataArrayIdType *revDescIndx, DimM1DescNbrer nbrer) const
172 if(!desc || !descIndx || !revDesc || !revDescIndx)
173 throw INTERP_KERNEL::Exception("MEDCouplingUMesh::buildDescendingConnectivityGen : present of a null pointer in input !");
174 checkConnectivityFullyDefined();
175 mcIdType nbOfCells=getNumberOfCells();
176 mcIdType nbOfNodes=getNumberOfNodes();
177 MCAuto<DataArrayIdType> revNodalIndx=DataArrayIdType::New(); revNodalIndx->alloc(nbOfNodes+1,1); revNodalIndx->fillWithZero();
178 mcIdType *revNodalIndxPtr=revNodalIndx->getPointer();
179 const mcIdType *conn=_nodal_connec->getConstPointer();
180 const mcIdType *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<DataArrayIdType> revDesc2(DataArrayIdType::New()); revDesc2->reserve(2*nbOfCells);
187 mcIdType *descIndxPtr=descIndx->getPointer(); *descIndxPtr++=0;
188 for(mcIdType eltId=0;eltId<nbOfCells;eltId++,descIndxPtr++)
190 mcIdType pos=connIndex[eltId];
191 mcIdType 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<mcIdType> tmp=new mcIdType[posP1-pos];
196 for(unsigned i=0;i<nbOfSons;i++)
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++)
202 revNodalIndxPtr[tmp[k]+1]++;
203 ret->insertNextCell(cmsId,nbOfNodesSon,tmp);
204 revDesc2->pushBackSilent(eltId);
206 descIndxPtr[0]=descIndxPtr[-1]+ToIdType(nbOfSons);
208 mcIdType nbOfCellsM1=ret->getNumberOfCells();
209 std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<mcIdType>());
210 MCAuto<DataArrayIdType> revNodal=DataArrayIdType::New(); revNodal->alloc(revNodalIndx->back(),1);
211 std::fill(revNodal->getPointer(),revNodal->getPointer()+revNodalIndx->back(),-1);
212 mcIdType *revNodalPtr=revNodal->getPointer();
213 const mcIdType *connM1=ret->getNodalConnectivity()->getConstPointer();
214 const mcIdType *connIndexM1=ret->getNodalConnectivityIndex()->getConstPointer();
215 for(mcIdType eltId=0;eltId<nbOfCellsM1;eltId++)
217 const mcIdType *strtNdlConnOfCurCell=connM1+connIndexM1[eltId]+1;
218 const mcIdType *endNdlConnOfCurCell=connM1+connIndexM1[eltId+1];
219 for(const mcIdType *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<mcIdType>(),-1))=eltId;
224 DataArrayIdType *commonCells=0,*commonCellsI=0;
225 FindCommonCellsAlg(3,0,ret->getNodalConnectivity(),ret->getNodalConnectivityIndex(),revNodal,revNodalIndx,commonCells,commonCellsI);
226 MCAuto<DataArrayIdType> commonCellsTmp(commonCells),commonCellsITmp(commonCellsI);
227 const mcIdType *commonCellsPtr(commonCells->getConstPointer()),*commonCellsIPtr(commonCellsI->getConstPointer());
228 mcIdType newNbOfCellsM1=-1;
229 MCAuto<DataArrayIdType> o2nM1=DataArrayIdType::ConvertIndexArrayToO2N(nbOfCellsM1,commonCells->begin(),
230 commonCellsI->begin(),commonCellsI->end(),newNbOfCellsM1);
231 std::vector<bool> isImpacted(nbOfCellsM1,false);
232 for(const mcIdType *work=commonCellsI->begin();work!=commonCellsI->end()-1;work++)
233 for(mcIdType work2=work[0];work2!=work[1];work2++)
234 isImpacted[commonCellsPtr[work2]]=true;
235 const mcIdType *o2nM1Ptr=o2nM1->getConstPointer();
236 MCAuto<DataArrayIdType> n2oM1=o2nM1->invertArrayO2N2N2OBis(newNbOfCellsM1);
237 const mcIdType *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 mcIdType *descPtr=desc->getPointer();
242 const INTERP_KERNEL::CellModel& cmsDft=INTERP_KERNEL::CellModel::GetCellModel(INTERP_KERNEL::NORM_POINT1);
243 for(mcIdType i=0;i<nbOfCellsM1;i++,descPtr++)
246 *descPtr=nbrer(o2nM1Ptr[i],0,cmsDft,false,0,0);
249 if(i!=n2oM1Ptr[o2nM1Ptr[i]])
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);
255 *descPtr=nbrer(o2nM1Ptr[i],0,cmsDft,false,0,0);
258 revDesc->reserve(newNbOfCellsM1);
259 revDescIndx->alloc(newNbOfCellsM1+1,1);
260 mcIdType *revDescIndxPtr=revDescIndx->getPointer(); *revDescIndxPtr++=0;
261 const mcIdType *revDesc2Ptr=revDesc2->getConstPointer();
262 for(mcIdType i=0;i<newNbOfCellsM1;i++,revDescIndxPtr++)
264 mcIdType oldCellIdM1=n2oM1Ptr[i];
265 if(!isImpacted[oldCellIdM1])
267 revDesc->pushBackSilent(revDesc2Ptr[oldCellIdM1]);
268 revDescIndxPtr[0]=revDescIndxPtr[-1]+1;
272 for(mcIdType j=commonCellsIPtr[0];j<commonCellsIPtr[1];j++)
273 revDesc->pushBackSilent(revDesc2Ptr[commonCellsPtr[j]]);
274 revDescIndxPtr[0]=revDescIndxPtr[-1]+commonCellsIPtr[1]-commonCellsIPtr[0];