1 // Copyright (C) 2007-2013 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.
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 (CEA/DEN)
21 #include "MEDCouplingExtrudedMesh.hxx"
22 #include "MEDCouplingUMesh.hxx"
23 #include "MEDCouplingMemArray.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
25 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
26 #include "CellModel.hxx"
28 #include "InterpolationUtils.hxx"
39 using namespace ParaMEDMEM;
42 * Build an extruded mesh instance from 3D and 2D unstructured mesh lying on the \b same \b coords.
43 * @param mesh3D 3D unstructured mesh.
44 * @param mesh2D 2D unstructured mesh lying on the same coordinates than mesh3D. \b Warning mesh2D is \b not \b const
45 * because the mesh is aggregated and potentially modified by rotate or translate method.
46 * @param cell2DId Id of cell in mesh2D mesh where the computation of 1D mesh will be done.
48 MEDCouplingExtrudedMesh *MEDCouplingExtrudedMesh::New(const MEDCouplingUMesh *mesh3D, const MEDCouplingUMesh *mesh2D, int cell2DId) throw(INTERP_KERNEL::Exception)
50 return new MEDCouplingExtrudedMesh(mesh3D,mesh2D,cell2DId);
54 * This constructor is here only for unserialisation process.
55 * This constructor is normally completely useless for end user.
57 MEDCouplingExtrudedMesh *MEDCouplingExtrudedMesh::New()
59 return new MEDCouplingExtrudedMesh;
62 MEDCouplingMeshType MEDCouplingExtrudedMesh::getType() const
67 std::size_t MEDCouplingExtrudedMesh::getHeapMemorySize() const
71 ret+=_mesh2D->getHeapMemorySize();
73 ret+=_mesh1D->getHeapMemorySize();
75 ret+=_mesh3D_ids->getHeapMemorySize();
76 return MEDCouplingMesh::getHeapMemorySize()+ret;
80 * This method copyies all tiny strings from other (name and components name).
81 * @throw if other and this have not same mesh type.
83 void MEDCouplingExtrudedMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
85 const MEDCouplingExtrudedMesh *otherC=dynamic_cast<const MEDCouplingExtrudedMesh *>(other);
87 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::copyTinyStringsFrom : meshes have not same type !");
88 MEDCouplingMesh::copyTinyStringsFrom(other);
89 _mesh2D->copyTinyStringsFrom(otherC->_mesh2D);
90 _mesh1D->copyTinyStringsFrom(otherC->_mesh1D);
93 MEDCouplingExtrudedMesh::MEDCouplingExtrudedMesh(const MEDCouplingUMesh *mesh3D, const MEDCouplingUMesh *mesh2D, int cell2DId) throw(INTERP_KERNEL::Exception)
94 try:_mesh2D(const_cast<MEDCouplingUMesh *>(mesh2D)),_mesh1D(MEDCouplingUMesh::New()),_mesh3D_ids(0),_cell_2D_id(cell2DId)
98 computeExtrusion(mesh3D);
99 setName(mesh3D->getName());
101 catch(INTERP_KERNEL::Exception& e)
108 _mesh3D_ids->decrRef();
112 MEDCouplingExtrudedMesh::MEDCouplingExtrudedMesh():_mesh2D(0),_mesh1D(0),_mesh3D_ids(0),_cell_2D_id(-1)
116 MEDCouplingExtrudedMesh::MEDCouplingExtrudedMesh(const MEDCouplingExtrudedMesh& other, bool deepCopy):MEDCouplingMesh(other),_cell_2D_id(other._cell_2D_id)
120 _mesh2D=other._mesh2D->clone(true);
121 _mesh1D=other._mesh1D->clone(true);
122 _mesh3D_ids=other._mesh3D_ids->deepCpy();
126 _mesh2D=other._mesh2D;
129 _mesh1D=other._mesh1D;
132 _mesh3D_ids=other._mesh3D_ids;
134 _mesh3D_ids->incrRef();
138 int MEDCouplingExtrudedMesh::getNumberOfCells() const
140 return _mesh2D->getNumberOfCells()*_mesh1D->getNumberOfCells();
143 int MEDCouplingExtrudedMesh::getNumberOfNodes() const
145 return _mesh2D->getNumberOfNodes();
148 int MEDCouplingExtrudedMesh::getSpaceDimension() const
153 int MEDCouplingExtrudedMesh::getMeshDimension() const
158 MEDCouplingMesh *MEDCouplingExtrudedMesh::deepCpy() const
163 MEDCouplingExtrudedMesh *MEDCouplingExtrudedMesh::clone(bool recDeepCpy) const
165 return new MEDCouplingExtrudedMesh(*this,recDeepCpy);
168 bool MEDCouplingExtrudedMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
171 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::isEqualIfNotWhy : input other pointer is null !");
172 const MEDCouplingExtrudedMesh *otherC=dynamic_cast<const MEDCouplingExtrudedMesh *>(other);
173 std::ostringstream oss;
176 reason="mesh given in input is not castable in MEDCouplingExtrudedMesh !";
179 if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
181 if(!_mesh2D->isEqualIfNotWhy(otherC->_mesh2D,prec,reason))
183 reason.insert(0,"Mesh2D unstructured meshes differ : ");
186 if(!_mesh1D->isEqualIfNotWhy(otherC->_mesh1D,prec,reason))
188 reason.insert(0,"Mesh1D unstructured meshes differ : ");
191 if(!_mesh3D_ids->isEqualIfNotWhy(*otherC->_mesh3D_ids,reason))
193 reason.insert(0,"Mesh3D ids DataArrayInt instances differ : ");
196 if(_cell_2D_id!=otherC->_cell_2D_id)
198 oss << "Cell 2D id of the two extruded mesh differ : this = " << _cell_2D_id << " other = " << otherC->_cell_2D_id;
205 bool MEDCouplingExtrudedMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
207 const MEDCouplingExtrudedMesh *otherC=dynamic_cast<const MEDCouplingExtrudedMesh *>(other);
210 if(!_mesh2D->isEqualWithoutConsideringStr(otherC->_mesh2D,prec))
212 if(!_mesh1D->isEqualWithoutConsideringStr(otherC->_mesh1D,prec))
214 if(!_mesh3D_ids->isEqualWithoutConsideringStr(*otherC->_mesh3D_ids))
216 if(_cell_2D_id!=otherC->_cell_2D_id)
221 void MEDCouplingExtrudedMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
222 DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
224 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::checkDeepEquivalWith : not implemented yet !");
227 void MEDCouplingExtrudedMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
228 DataArrayInt *&cellCor) const throw(INTERP_KERNEL::Exception)
230 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::checkDeepEquivalOnSameNodesWith : not implemented yet !");
233 INTERP_KERNEL::NormalizedCellType MEDCouplingExtrudedMesh::getTypeOfCell(int cellId) const
235 const int *ids=_mesh3D_ids->getConstPointer();
236 int nbOf3DCells=_mesh3D_ids->getNumberOfTuples();
237 const int *where=std::find(ids,ids+nbOf3DCells,cellId);
238 if(where==ids+nbOf3DCells)
239 throw INTERP_KERNEL::Exception("Invalid cellId specified >= getNumberOfCells() !");
240 int nbOfCells2D=_mesh2D->getNumberOfCells();
241 int locId=((int)std::distance(ids,where))%nbOfCells2D;
242 INTERP_KERNEL::NormalizedCellType tmp=_mesh2D->getTypeOfCell(locId);
243 return INTERP_KERNEL::CellModel::GetCellModel(tmp).getExtrudedType();
246 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingExtrudedMesh::getAllGeoTypes() const
248 const std::set<INTERP_KERNEL::NormalizedCellType>& ret2D=_mesh2D->getAllTypes();
249 std::set<INTERP_KERNEL::NormalizedCellType> ret;
250 for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=ret2D.begin();it!=ret2D.end();it++)
251 ret.insert(INTERP_KERNEL::CellModel::GetCellModel(*it).getExtrudedType());
255 DataArrayInt *MEDCouplingExtrudedMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
257 const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
258 INTERP_KERNEL::NormalizedCellType revExtTyp=cm.getReverseExtrudedType();
259 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
260 if(revExtTyp==INTERP_KERNEL::NORM_ERROR)
265 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp=_mesh2D->giveCellsWithType(revExtTyp);
266 int nbOfLevs=_mesh1D->getNumberOfCells();
267 int nbOfCells2D=_mesh2D->getNumberOfCells();
268 int nbOfTuples=tmp->getNumberOfTuples();
269 ret->alloc(nbOfLevs*nbOfTuples,1);
270 int *pt=ret->getPointer();
271 for(int i=0;i<nbOfLevs;i++,pt+=nbOfTuples)
272 std::transform(tmp->begin(),tmp->end(),pt,std::bind2nd(std::plus<double>(),i*nbOfCells2D));
273 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2=ret->renumberR(_mesh3D_ids->begin());
278 DataArrayInt *MEDCouplingExtrudedMesh::computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
280 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2D=_mesh2D->computeNbOfNodesPerCell();
281 int nbOfLevs=_mesh1D->getNumberOfCells();
282 int nbOfCells2D=_mesh2D->getNumberOfCells();
283 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret3D=DataArrayInt::New(); ret3D->alloc(nbOfLevs*nbOfCells2D,1);
284 int *pt=ret3D->getPointer();
285 for(int i=0;i<nbOfLevs;i++,pt+=nbOfCells2D)
286 std::copy(ret2D->begin(),ret2D->end(),pt);
287 ret3D->applyLin(2,0,0);
288 return ret3D->renumberR(_mesh3D_ids->begin());
291 DataArrayInt *MEDCouplingExtrudedMesh::computeNbOfFacesPerCell() const throw(INTERP_KERNEL::Exception)
293 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2D=_mesh2D->computeNbOfNodesPerCell();
294 int nbOfLevs=_mesh1D->getNumberOfCells();
295 int nbOfCells2D=_mesh2D->getNumberOfCells();
296 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret3D=DataArrayInt::New(); ret3D->alloc(nbOfLevs*nbOfCells2D,1);
297 int *pt=ret3D->getPointer();
298 for(int i=0;i<nbOfLevs;i++,pt+=nbOfCells2D)
299 std::copy(ret2D->begin(),ret2D->end(),pt);
300 ret3D->applyLin(2,2,0);
301 return ret3D->renumberR(_mesh3D_ids->begin());
304 int MEDCouplingExtrudedMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
307 int nbOfCells2D=_mesh2D->getNumberOfCells();
308 for(int i=0;i<nbOfCells2D;i++)
310 INTERP_KERNEL::NormalizedCellType t=_mesh2D->getTypeOfCell(i);
311 if(INTERP_KERNEL::CellModel::GetCellModel(t).getExtrudedType()==type)
314 return ret*_mesh1D->getNumberOfCells();
317 void MEDCouplingExtrudedMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
319 int nbOfCells2D=_mesh2D->getNumberOfCells();
320 int nbOfNodes2D=_mesh2D->getNumberOfNodes();
321 int locId=cellId%nbOfCells2D;
322 int lev=cellId/nbOfCells2D;
323 std::vector<int> tmp,tmp2;
324 _mesh2D->getNodeIdsOfCell(locId,tmp);
326 std::transform(tmp.begin(),tmp.end(),tmp.begin(),std::bind2nd(std::plus<int>(),nbOfNodes2D*lev));
327 std::transform(tmp2.begin(),tmp2.end(),tmp2.begin(),std::bind2nd(std::plus<int>(),nbOfNodes2D*(lev+1)));
328 conn.insert(conn.end(),tmp.begin(),tmp.end());
329 conn.insert(conn.end(),tmp2.begin(),tmp2.end());
332 void MEDCouplingExtrudedMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
334 int nbOfNodes2D=_mesh2D->getNumberOfNodes();
335 int locId=nodeId%nbOfNodes2D;
336 int lev=nodeId/nbOfNodes2D;
337 std::vector<double> tmp,tmp2;
338 _mesh2D->getCoordinatesOfNode(locId,tmp);
340 int spaceDim=_mesh1D->getSpaceDimension();
341 const double *z=_mesh1D->getCoords()->getConstPointer();
342 std::transform(tmp.begin(),tmp.end(),z+lev*spaceDim,tmp.begin(),std::plus<double>());
343 std::transform(tmp2.begin(),tmp2.end(),z+(lev+1)*spaceDim,tmp2.begin(),std::plus<double>());
344 coo.insert(coo.end(),tmp.begin(),tmp.end());
345 coo.insert(coo.end(),tmp2.begin(),tmp2.end());
348 std::string MEDCouplingExtrudedMesh::simpleRepr() const
350 std::ostringstream ret;
351 ret << "3D Extruded mesh from a 2D Surf Mesh with name : \"" << getName() << "\"\n";
352 ret << "Description of mesh : \"" << getDescription() << "\"\n";
354 double tt=getTime(tmpp1,tmpp2);
355 ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
356 ret << "Iteration : " << tmpp1 << " Order : " << tmpp2 << "\n";
357 ret << "Cell id where 1D mesh has been deduced : " << _cell_2D_id << "\n";
358 ret << "Number of cells : " << getNumberOfCells() << "(" << _mesh2D->getNumberOfCells() << "x" << _mesh1D->getNumberOfCells() << ")\n";
359 ret << "1D Mesh info : _____________________\n\n\n";
360 ret << _mesh1D->simpleRepr();
361 ret << "\n\n\n2D Mesh info : _____________________\n\n\n" << _mesh2D->simpleRepr() << "\n\n\n";
365 std::string MEDCouplingExtrudedMesh::advancedRepr() const
367 std::ostringstream ret;
368 ret << "3D Extruded mesh from a 2D Surf Mesh with name : \"" << getName() << "\"\n";
369 ret << "Description of mesh : \"" << getDescription() << "\"\n";
371 double tt=getTime(tmpp1,tmpp2);
372 ret << "Time attached to the mesh (unit) : " << tt << " (" << getTimeUnit() << ")\n";
373 ret << "Iteration : " << tmpp1 << " Order : " << tmpp2 << "\n";
374 ret << "Cell id where 1D mesh has been deduced : " << _cell_2D_id << "\n";
375 ret << "Number of cells : " << getNumberOfCells() << "(" << _mesh2D->getNumberOfCells() << "x" << _mesh1D->getNumberOfCells() << ")\n";
376 ret << "1D Mesh info : _____________________\n\n\n";
377 ret << _mesh1D->advancedRepr();
378 ret << "\n\n\n2D Mesh info : _____________________\n\n\n" << _mesh2D->advancedRepr() << "\n\n\n";
379 ret << "3D cell ids per level :\n";
383 void MEDCouplingExtrudedMesh::checkCoherency() const throw (INTERP_KERNEL::Exception)
387 void MEDCouplingExtrudedMesh::checkCoherency1(double eps) const throw(INTERP_KERNEL::Exception)
392 void MEDCouplingExtrudedMesh::checkCoherency2(double eps) const throw(INTERP_KERNEL::Exception)
394 checkCoherency1(eps);
397 void MEDCouplingExtrudedMesh::getBoundingBox(double *bbox) const
400 _mesh2D->getBoundingBox(bbox2D);
401 const double *nodes1D=_mesh1D->getCoords()->getConstPointer();
402 int nbOfNodes1D=_mesh1D->getNumberOfNodes();
403 double bbox1DMin[3],bbox1DMax[3],tmp[3];
404 std::fill(bbox1DMin,bbox1DMin+3,std::numeric_limits<double>::max());
405 std::fill(bbox1DMax,bbox1DMax+3,-(std::numeric_limits<double>::max()));
406 for(int i=0;i<nbOfNodes1D;i++)
408 std::transform(nodes1D+3*i,nodes1D+3*(i+1),bbox1DMin,bbox1DMin,static_cast<const double& (*)(const double&, const double&)>(std::min<double>));
409 std::transform(nodes1D+3*i,nodes1D+3*(i+1),bbox1DMax,bbox1DMax,static_cast<const double& (*)(const double&, const double&)>(std::max<double>));
411 std::transform(bbox1DMax,bbox1DMax+3,bbox1DMin,tmp,std::minus<double>());
412 int id=(int)std::distance(tmp,std::max_element(tmp,tmp+3));
413 bbox[0]=bbox1DMin[0]; bbox[1]=bbox1DMax[0];
414 bbox[2]=bbox1DMin[1]; bbox[3]=bbox1DMax[1];
415 bbox[4]=bbox1DMin[2]; bbox[5]=bbox1DMax[2];
416 bbox[2*id+1]+=tmp[id];
419 void MEDCouplingExtrudedMesh::updateTime() const
423 updateTimeWith(*_mesh2D);
427 updateTimeWith(*_mesh1D);
431 void MEDCouplingExtrudedMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
433 throw INTERP_KERNEL::Exception("Functionnality of renumbering cells unavailable for ExtrudedMesh");
436 MEDCouplingUMesh *MEDCouplingExtrudedMesh::build3DUnstructuredMesh() const
438 MEDCouplingUMesh *ret=_mesh2D->buildExtrudedMesh(_mesh1D,0);
439 const int *renum=_mesh3D_ids->getConstPointer();
440 ret->renumberCells(renum,false);
441 ret->setName(getName());
445 MEDCouplingUMesh *MEDCouplingExtrudedMesh::buildUnstructured() const throw(INTERP_KERNEL::Exception)
447 return build3DUnstructuredMesh();
450 MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::getMeasureField(bool) const
452 std::string name="MeasureOfMesh_";
454 MEDCouplingFieldDouble *ret2D=_mesh2D->getMeasureField(true);
455 MEDCouplingFieldDouble *ret1D=_mesh1D->getMeasureField(true);
456 const double *ret2DPtr=ret2D->getArray()->getConstPointer();
457 const double *ret1DPtr=ret1D->getArray()->getConstPointer();
458 int nbOf2DCells=_mesh2D->getNumberOfCells();
459 int nbOf1DCells=_mesh1D->getNumberOfCells();
460 int nbOf3DCells=nbOf2DCells*nbOf1DCells;
461 const int *renum=_mesh3D_ids->getConstPointer();
462 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
464 ret->synchronizeTimeWithMesh();
465 DataArrayDouble *da=DataArrayDouble::New();
466 da->alloc(nbOf3DCells,1);
467 double *retPtr=da->getPointer();
468 for(int i=0;i<nbOf1DCells;i++)
469 for(int j=0;j<nbOf2DCells;j++)
470 retPtr[renum[i*nbOf2DCells+j]]=ret2DPtr[j]*ret1DPtr[i];
473 ret->setName(name.c_str());
479 MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::getMeasureFieldOnNode(bool isAbs) const
481 //not implemented yet
485 MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::buildOrthogonalField() const
487 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::buildOrthogonalField : This method has no sense for MEDCouplingExtrudedMesh that is 3D !");
490 int MEDCouplingExtrudedMesh::getCellContainingPoint(const double *pos, double eps) const
492 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::getCellContainingPoint : not implemented yet !");
495 MEDCouplingExtrudedMesh::~MEDCouplingExtrudedMesh()
502 _mesh3D_ids->decrRef();
505 void MEDCouplingExtrudedMesh::computeExtrusion(const MEDCouplingUMesh *mesh3D) throw(INTERP_KERNEL::Exception)
507 const char errMsg1[]="2D mesh is empty unable to compute extrusion !";
508 const char errMsg2[]="Coords between 2D and 3D meshes are not the same ! Try MEDCouplingPointSet::tryToShareSameCoords method";
509 const char errMsg3[]="No chance to find extrusion pattern in mesh3D,mesh2D couple because nbCells3D%nbCells2D!=0 !";
510 if(_mesh2D==0 || mesh3D==0)
511 throw INTERP_KERNEL::Exception(errMsg1);
512 if(_mesh2D->getCoords()!=mesh3D->getCoords())
513 throw INTERP_KERNEL::Exception(errMsg2);
514 if(mesh3D->getNumberOfCells()%_mesh2D->getNumberOfCells()!=0)
515 throw INTERP_KERNEL::Exception(errMsg3);
517 _mesh3D_ids=DataArrayInt::New();
519 _mesh1D=MEDCouplingUMesh::New();
520 computeExtrusionAlg(mesh3D);
523 void MEDCouplingExtrudedMesh::build1DExtrusion(int idIn3DDesc, int newId, int nbOf1DLev, MEDCouplingUMesh *subMesh,
524 const int *desc3D, const int *descIndx3D,
525 const int *revDesc3D, const int *revDescIndx3D,
526 bool computeMesh1D) throw(INTERP_KERNEL::Exception)
528 int nbOf2DCells=_mesh2D->getNumberOfCells();
529 int start=revDescIndx3D[idIn3DDesc];
530 int end=revDescIndx3D[idIn3DDesc+1];
533 std::ostringstream ost; ost << "Invalid bases 2D mesh specified : 2D cell # " << idIn3DDesc;
534 ost << " shared by more than 1 3D cell !!!";
535 throw INTERP_KERNEL::Exception(ost.str().c_str());
537 int current3DCell=revDesc3D[start];
538 int current2DCell=idIn3DDesc;
539 int *mesh3DIDs=_mesh3D_ids->getPointer();
540 mesh3DIDs[newId]=current3DCell;
541 const int *conn2D=subMesh->getNodalConnectivity()->getConstPointer();
542 const int *conn2DIndx=subMesh->getNodalConnectivityIndex()->getConstPointer();
543 for(int i=1;i<nbOf1DLev;i++)
545 std::vector<int> conn(conn2D+conn2DIndx[current2DCell]+1,conn2D+conn2DIndx[current2DCell+1]);
546 std::sort(conn.begin(),conn.end());
548 computeBaryCenterOfFace(conn,i-1);
549 current2DCell=findOppositeFaceOf(current2DCell,current3DCell,conn,
550 desc3D,descIndx3D,conn2D,conn2DIndx);
551 start=revDescIndx3D[current2DCell];
552 end=revDescIndx3D[current2DCell+1];
555 std::ostringstream ost; ost << "Expecting to have 2 3D cells attached to 2D cell " << current2DCell << "!";
556 ost << " : Impossible or call tryToShareSameCoords method !";
557 throw INTERP_KERNEL::Exception(ost.str().c_str());
559 if(revDesc3D[start]!=current3DCell)
560 current3DCell=revDesc3D[start];
562 current3DCell=revDesc3D[start+1];
563 mesh3DIDs[i*nbOf2DCells+newId]=current3DCell;
567 std::vector<int> conn(conn2D+conn2DIndx[current2DCell]+1,conn2D+conn2DIndx[current2DCell+1]);
568 std::sort(conn.begin(),conn.end());
569 computeBaryCenterOfFace(conn,nbOf1DLev-1);
570 current2DCell=findOppositeFaceOf(current2DCell,current3DCell,conn,
571 desc3D,descIndx3D,conn2D,conn2DIndx);
573 conn.insert(conn.end(),conn2D+conn2DIndx[current2DCell]+1,conn2D+conn2DIndx[current2DCell+1]);
574 std::sort(conn.begin(),conn.end());
575 computeBaryCenterOfFace(conn,nbOf1DLev);
579 int MEDCouplingExtrudedMesh::findOppositeFaceOf(int current2DCell, int current3DCell, const std::vector<int>& connSorted,
580 const int *desc3D, const int *descIndx3D,
581 const int *conn2D, const int *conn2DIndx) throw(INTERP_KERNEL::Exception)
583 int start=descIndx3D[current3DCell];
584 int end=descIndx3D[current3DCell+1];
586 for(const int *candidate2D=desc3D+start;candidate2D!=desc3D+end && !found;candidate2D++)
588 if(*candidate2D!=current2DCell)
590 std::vector<int> conn2(conn2D+conn2DIndx[*candidate2D]+1,conn2D+conn2DIndx[*candidate2D+1]);
591 std::sort(conn2.begin(),conn2.end());
592 std::list<int> intersect;
593 std::set_intersection(connSorted.begin(),connSorted.end(),conn2.begin(),conn2.end(),
594 std::insert_iterator< std::list<int> >(intersect,intersect.begin()));
595 if(intersect.empty())
599 std::ostringstream ost; ost << "Impossible to find an opposite 2D face of face # " << current2DCell;
600 ost << " in 3D cell # " << current3DCell << " : Impossible or call tryToShareSameCoords method !";
601 throw INTERP_KERNEL::Exception(ost.str().c_str());
604 void MEDCouplingExtrudedMesh::computeBaryCenterOfFace(const std::vector<int>& nodalConnec, int lev1DId)
606 double *zoneToUpdate=_mesh1D->getCoords()->getPointer()+lev1DId*3;
607 std::fill(zoneToUpdate,zoneToUpdate+3,0.);
608 const double *coords=_mesh2D->getCoords()->getConstPointer();
609 for(std::vector<int>::const_iterator iter=nodalConnec.begin();iter!=nodalConnec.end();iter++)
610 std::transform(zoneToUpdate,zoneToUpdate+3,coords+3*(*iter),zoneToUpdate,std::plus<double>());
611 std::transform(zoneToUpdate,zoneToUpdate+3,zoneToUpdate,std::bind2nd(std::multiplies<double>(),(double)(1./(int)nodalConnec.size())));
614 int MEDCouplingExtrudedMesh::FindCorrespCellByNodalConn(const std::vector<int>& nodalConnec, const int *revNodalPtr, const int *revNodalIndxPtr) throw(INTERP_KERNEL::Exception)
616 std::vector<int>::const_iterator iter=nodalConnec.begin();
617 std::set<int> s1(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1]);
619 for(;iter!=nodalConnec.end();iter++)
621 std::set<int> s2(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1]);
623 std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(),std::insert_iterator< std::set<int> >(s3,s3.end()));
627 return *(s1.begin());
628 std::ostringstream ostr;
629 ostr << "Cell with nodal connec : ";
630 std::copy(nodalConnec.begin(),nodalConnec.end(),std::ostream_iterator<int>(ostr," "));
631 ostr << " is not part of mesh";
632 throw INTERP_KERNEL::Exception(ostr.str().c_str());
636 * This method is callable on 1Dmeshes (meshDim==1 && spaceDim==3) returned by MEDCouplingExtrudedMesh::getMesh1D typically.
637 * These 1Dmeshes (meshDim==1 && spaceDim==3) have a special semantic because these meshes do not specify a static location but a translation along a path.
638 * This method checks that 'm1' and 'm2' are compatible, if not an exception is thrown. In case these meshes ('m1' and 'm2') are compatible 2 corresponding meshes
639 * are created ('m1r' and 'm2r') that can be used for interpolation.
640 * @param m1 input mesh with meshDim==1 and spaceDim==3
641 * @param m2 input mesh with meshDim==1 and spaceDim==3
642 * @param eps tolerance acceptable to determine compatibility
643 * @param m1r output mesh with ref count equal to 1 with meshDim==1 and spaceDim==1
644 * @param m2r output mesh with ref count equal to 1 with meshDim==1 and spaceDim==1
645 * @param v is the output normalized vector of the common direction of 'm1' and 'm2'
646 * @throw in case that m1 and m2 are not compatible each other.
648 void MEDCouplingExtrudedMesh::Project1DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
649 MEDCouplingUMesh *&m1r, MEDCouplingUMesh *&m2r, double *v) throw(INTERP_KERNEL::Exception)
651 if(m1->getSpaceDimension()!=3 || m1->getSpaceDimension()!=3)
652 throw INTERP_KERNEL::Exception("Input meshes are expected to have a spaceDim==3 for Projec1D !");
655 m1r->changeSpaceDimension(1);
656 m2r->changeSpaceDimension(1);
658 std::vector<double> ref,ref2;
659 m1->getNodeIdsOfCell(0,c);
660 m1->getCoordinatesOfNode(c[0],ref);
661 m1->getCoordinatesOfNode(c[1],ref2);
662 std::transform(ref2.begin(),ref2.end(),ref.begin(),v,std::minus<double>());
663 double n=INTERP_KERNEL::norm<3>(v);
664 std::transform(v,v+3,v,std::bind2nd(std::multiplies<double>(),1/n));
665 m1->project1D(&ref[0],v,eps,m1r->getCoords()->getPointer());
666 m2->project1D(&ref[0],v,eps,m2r->getCoords()->getPointer());
670 void MEDCouplingExtrudedMesh::rotate(const double *center, const double *vector, double angle)
672 _mesh2D->rotate(center,vector,angle);
673 _mesh1D->rotate(center,vector,angle);
676 void MEDCouplingExtrudedMesh::translate(const double *vector)
678 _mesh2D->translate(vector);
679 _mesh1D->translate(vector);
682 void MEDCouplingExtrudedMesh::scale(const double *point, double factor)
684 _mesh2D->scale(point,factor);
685 _mesh1D->scale(point,factor);
688 std::vector<int> MEDCouplingExtrudedMesh::getDistributionOfTypes() const throw(INTERP_KERNEL::Exception)
690 throw INTERP_KERNEL::Exception("Not implemented yet !");
693 DataArrayInt *MEDCouplingExtrudedMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
695 throw INTERP_KERNEL::Exception("Not implemented yet !");
698 void MEDCouplingExtrudedMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
700 throw INTERP_KERNEL::Exception("Not implemented yet !");
703 MEDCouplingMesh *MEDCouplingExtrudedMesh::buildPart(const int *start, const int *end) const
705 // not implemented yet !
709 MEDCouplingMesh *MEDCouplingExtrudedMesh::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
711 // not implemented yet !
715 DataArrayInt *MEDCouplingExtrudedMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception)
717 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::simplexize : unavailable for such a type of mesh : Extruded !");
720 MEDCouplingMesh *MEDCouplingExtrudedMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
722 // not implemented yet !
726 DataArrayDouble *MEDCouplingExtrudedMesh::getCoordinatesAndOwner() const
728 DataArrayDouble *arr2D=_mesh2D->getCoords();
729 DataArrayDouble *arr1D=_mesh1D->getCoords();
730 DataArrayDouble *ret=DataArrayDouble::New();
731 ret->alloc(getNumberOfNodes(),3);
732 int nbOf1DLev=_mesh1D->getNumberOfNodes();
733 int nbOf2DNodes=_mesh2D->getNumberOfNodes();
734 const double *ptSrc=arr2D->getConstPointer();
735 double *pt=ret->getPointer();
736 std::copy(ptSrc,ptSrc+3*nbOf2DNodes,pt);
737 for(int i=1;i<nbOf1DLev;i++)
739 std::copy(ptSrc,ptSrc+3*nbOf2DNodes,pt+3*i*nbOf2DNodes);
741 std::copy(arr1D->getConstPointer()+3*i,arr1D->getConstPointer()+3*(i+1),vec);
742 std::transform(arr1D->getConstPointer()+3*(i-1),arr1D->getConstPointer()+3*i,vec,vec,std::minus<double>());
743 for(int j=0;j<nbOf2DNodes;j++)
744 std::transform(vec,vec+3,pt+3*(i*nbOf2DNodes+j),pt+3*(i*nbOf2DNodes+j),std::plus<double>());
749 DataArrayDouble *MEDCouplingExtrudedMesh::getBarycenterAndOwner() const
751 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::getBarycenterAndOwner : not yet implemented !");
754 DataArrayDouble *MEDCouplingExtrudedMesh::computeIsoBarycenterOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
756 throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::computeIsoBarycenterOfNodesPerCell: not yet implemented !");
759 void MEDCouplingExtrudedMesh::computeExtrusionAlg(const MEDCouplingUMesh *mesh3D) throw(INTERP_KERNEL::Exception)
761 _mesh3D_ids->alloc(mesh3D->getNumberOfCells(),1);
762 int nbOf1DLev=mesh3D->getNumberOfCells()/_mesh2D->getNumberOfCells();
763 _mesh1D->setMeshDimension(1);
764 _mesh1D->allocateCells(nbOf1DLev);
766 for(int i=0;i<nbOf1DLev;i++)
770 _mesh1D->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,tmpConn);
772 _mesh1D->finishInsertingCells();
773 DataArrayDouble *myCoords=DataArrayDouble::New();
774 myCoords->alloc(nbOf1DLev+1,3);
775 _mesh1D->setCoords(myCoords);
777 DataArrayInt *desc,*descIndx,*revDesc,*revDescIndx;
778 desc=DataArrayInt::New(); descIndx=DataArrayInt::New(); revDesc=DataArrayInt::New(); revDescIndx=DataArrayInt::New();
779 MEDCouplingUMesh *subMesh=mesh3D->buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
780 DataArrayInt *revNodal2D,*revNodalIndx2D;
781 revNodal2D=DataArrayInt::New(); revNodalIndx2D=DataArrayInt::New();
782 subMesh->getReverseNodalConnectivity(revNodal2D,revNodalIndx2D);
783 const int *nodal2D=_mesh2D->getNodalConnectivity()->getConstPointer();
784 const int *nodal2DIndx=_mesh2D->getNodalConnectivityIndex()->getConstPointer();
785 const int *revNodal2DPtr=revNodal2D->getConstPointer();
786 const int *revNodalIndx2DPtr=revNodalIndx2D->getConstPointer();
787 const int *descP=desc->getConstPointer();
788 const int *descIndxP=descIndx->getConstPointer();
789 const int *revDescP=revDesc->getConstPointer();
790 const int *revDescIndxP=revDescIndx->getConstPointer();
792 int nbOf2DCells=_mesh2D->getNumberOfCells();
793 for(int i=0;i<nbOf2DCells;i++)
796 std::vector<int> nodalConnec(nodal2D+nodal2DIndx[i]+1,nodal2D+nodal2DIndx[i+1]);
799 idInSubMesh=FindCorrespCellByNodalConn(nodalConnec,revNodal2DPtr,revNodalIndx2DPtr);
801 catch(INTERP_KERNEL::Exception& e)
803 std::ostringstream ostr; ostr << "mesh2D cell # " << i << " is not part of any cell of 3D mesh !\n";
805 throw INTERP_KERNEL::Exception(ostr.str().c_str());
807 build1DExtrusion(idInSubMesh,i,nbOf1DLev,subMesh,descP,descIndxP,revDescP,revDescIndxP,i==_cell_2D_id);
810 revNodal2D->decrRef();
811 revNodalIndx2D->decrRef();
816 revDescIndx->decrRef();
819 void MEDCouplingExtrudedMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
821 std::vector<int> tinyInfo1;
822 std::vector<std::string> ls1;
823 std::vector<double> ls3;
824 _mesh2D->getTinySerializationInformation(ls3,tinyInfo1,ls1);
825 std::vector<int> tinyInfo2;
826 std::vector<std::string> ls2;
827 std::vector<double> ls4;
828 _mesh1D->getTinySerializationInformation(ls4,tinyInfo2,ls2);
829 tinyInfo.clear(); littleStrings.clear();
830 tinyInfo.insert(tinyInfo.end(),tinyInfo1.begin(),tinyInfo1.end());
831 littleStrings.insert(littleStrings.end(),ls1.begin(),ls1.end());
832 tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
833 littleStrings.insert(littleStrings.end(),ls2.begin(),ls2.end());
834 tinyInfo.push_back(_cell_2D_id);
835 tinyInfo.push_back((int)tinyInfo1.size());
836 tinyInfo.push_back(_mesh3D_ids->getNbOfElems());
837 littleStrings.push_back(getName());
838 littleStrings.push_back(getDescription());
841 void MEDCouplingExtrudedMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
843 std::size_t sz=tinyInfo.size();
844 int sz1=tinyInfo[sz-2];
845 std::vector<int> ti1(tinyInfo.begin(),tinyInfo.begin()+sz1);
846 std::vector<int> ti2(tinyInfo.begin()+sz1,tinyInfo.end()-3);
847 MEDCouplingUMesh *um=MEDCouplingUMesh::New();
848 DataArrayInt *a1tmp=DataArrayInt::New();
849 DataArrayDouble *a2tmp=DataArrayDouble::New();
851 std::vector<std::string> ls1,ls2;
852 um->resizeForUnserialization(ti1,a1tmp,a2tmp,ls1);
853 la1+=a1tmp->getNbOfElems(); la2+=a2tmp->getNbOfElems();
854 a1tmp->decrRef(); a2tmp->decrRef();
855 a1tmp=DataArrayInt::New(); a2tmp=DataArrayDouble::New();
856 um->resizeForUnserialization(ti2,a1tmp,a2tmp,ls2);
857 la1+=a1tmp->getNbOfElems(); la2+=a2tmp->getNbOfElems();
858 a1tmp->decrRef(); a2tmp->decrRef();
861 a1->alloc(la1+tinyInfo[sz-1],1);
863 littleStrings.resize(ls1.size()+ls2.size()+2);
866 void MEDCouplingExtrudedMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
868 a1=DataArrayInt::New(); a2=DataArrayDouble::New();
869 DataArrayInt *a1_1=0,*a1_2=0;
870 DataArrayDouble *a2_1=0,*a2_2=0;
871 _mesh2D->serialize(a1_1,a2_1);
872 _mesh1D->serialize(a1_2,a2_2);
873 a1->alloc(a1_1->getNbOfElems()+a1_2->getNbOfElems()+_mesh3D_ids->getNbOfElems(),1);
874 int *ptri=a1->getPointer();
875 ptri=std::copy(a1_1->getConstPointer(),a1_1->getConstPointer()+a1_1->getNbOfElems(),ptri);
877 ptri=std::copy(a1_2->getConstPointer(),a1_2->getConstPointer()+a1_2->getNbOfElems(),ptri);
879 std::copy(_mesh3D_ids->getConstPointer(),_mesh3D_ids->getConstPointer()+_mesh3D_ids->getNbOfElems(),ptri);
880 a2->alloc(a2_1->getNbOfElems()+a2_2->getNbOfElems(),1);
881 double *ptrd=a2->getPointer();
882 ptrd=std::copy(a2_1->getConstPointer(),a2_1->getConstPointer()+a2_1->getNbOfElems(),ptrd);
884 std::copy(a2_2->getConstPointer(),a2_2->getConstPointer()+a2_2->getNbOfElems(),ptrd);
888 void MEDCouplingExtrudedMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
890 setName(littleStrings[littleStrings.size()-2].c_str());
891 setDescription(littleStrings.back().c_str());
892 std::size_t sz=tinyInfo.size();
893 int sz1=tinyInfo[sz-2];
894 _cell_2D_id=tinyInfo[sz-3];
895 std::vector<int> ti1(tinyInfo.begin(),tinyInfo.begin()+sz1);
896 std::vector<int> ti2(tinyInfo.begin()+sz1,tinyInfo.end()-3);
897 DataArrayInt *a1tmp=DataArrayInt::New();
898 DataArrayDouble *a2tmp=DataArrayDouble::New();
899 const int *a1Ptr=a1->getConstPointer();
900 const double *a2Ptr=a2->getConstPointer();
901 _mesh2D=MEDCouplingUMesh::New();
902 std::vector<std::string> ls1,ls2;
903 _mesh2D->resizeForUnserialization(ti1,a1tmp,a2tmp,ls1);
904 std::copy(a2Ptr,a2Ptr+a2tmp->getNbOfElems(),a2tmp->getPointer());
905 std::copy(a1Ptr,a1Ptr+a1tmp->getNbOfElems(),a1tmp->getPointer());
906 a2Ptr+=a2tmp->getNbOfElems();
907 a1Ptr+=a1tmp->getNbOfElems();
908 ls2.insert(ls2.end(),littleStrings.begin(),littleStrings.begin()+ls1.size());
909 std::vector<double> d1(1);
910 _mesh2D->unserialization(d1,ti1,a1tmp,a2tmp,ls2);
911 a1tmp->decrRef(); a2tmp->decrRef();
914 ls2.insert(ls2.end(),littleStrings.begin()+ls1.size(),littleStrings.end()-2);
915 _mesh1D=MEDCouplingUMesh::New();
916 a1tmp=DataArrayInt::New(); a2tmp=DataArrayDouble::New();
917 _mesh1D->resizeForUnserialization(ti2,a1tmp,a2tmp,ls1);
918 std::copy(a2Ptr,a2Ptr+a2tmp->getNbOfElems(),a2tmp->getPointer());
919 std::copy(a1Ptr,a1Ptr+a1tmp->getNbOfElems(),a1tmp->getPointer());
920 a1Ptr+=a1tmp->getNbOfElems();
921 _mesh1D->unserialization(d1,ti2,a1tmp,a2tmp,ls2);
922 a1tmp->decrRef(); a2tmp->decrRef();
924 _mesh3D_ids=DataArrayInt::New();
925 int szIds=(int)std::distance(a1Ptr,a1->getConstPointer()+a1->getNbOfElems());
926 _mesh3D_ids->alloc(szIds,1);
927 std::copy(a1Ptr,a1Ptr+szIds,_mesh3D_ids->getPointer());
930 void MEDCouplingExtrudedMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception)
932 MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
933 m->writeVTKLL(ofs,cellData,pointData);
936 void MEDCouplingExtrudedMesh::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
938 stream << "MEDCouplingExtrudedMesh C++ instance at " << this << ". Name : \"" << getName() << "\".";
941 std::string MEDCouplingExtrudedMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception)
943 return _mesh2D->getVTKDataSetType();