Salome HOME
Merge from V6_main (04/10/2012)
[modules/med.git] / src / MEDCoupling / MEDCouplingExtrudedMesh.cxx
1 // Copyright (C) 2007-2012  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.
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 (CEA/DEN)
20
21 #include "MEDCouplingExtrudedMesh.hxx"
22 #include "MEDCouplingUMesh.hxx"
23 #include "MEDCouplingMemArray.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
25 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
26 #include "CellModel.hxx"
27
28 #include "InterpolationUtils.hxx"
29
30 #include <limits>
31 #include <algorithm>
32 #include <functional>
33 #include <iterator>
34 #include <sstream>
35 #include <cmath>
36 #include <list>
37 #include <set>
38
39 using namespace ParaMEDMEM;
40
41 /*!
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.
47  */
48 MEDCouplingExtrudedMesh *MEDCouplingExtrudedMesh::New(const MEDCouplingUMesh *mesh3D, const MEDCouplingUMesh *mesh2D, int cell2DId) throw(INTERP_KERNEL::Exception)
49 {
50   return new MEDCouplingExtrudedMesh(mesh3D,mesh2D,cell2DId);
51 }
52
53 /*!
54  * This constructor is here only for unserialisation process.
55  * This constructor is normally completely useless for end user.
56  */
57 MEDCouplingExtrudedMesh *MEDCouplingExtrudedMesh::New()
58 {
59   return new MEDCouplingExtrudedMesh;
60 }
61
62 MEDCouplingMeshType MEDCouplingExtrudedMesh::getType() const
63 {
64   return EXTRUDED;
65 }
66
67 /*!
68  * This method copyies all tiny strings from other (name and components name).
69  * @throw if other and this have not same mesh type.
70  */
71 void MEDCouplingExtrudedMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
72 {
73   const MEDCouplingExtrudedMesh *otherC=dynamic_cast<const MEDCouplingExtrudedMesh *>(other);
74   if(!otherC)
75     throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::copyTinyStringsFrom : meshes have not same type !");
76   MEDCouplingMesh::copyTinyStringsFrom(other);
77   _mesh2D->copyTinyStringsFrom(otherC->_mesh2D);
78   _mesh1D->copyTinyStringsFrom(otherC->_mesh1D);
79 }
80
81 MEDCouplingExtrudedMesh::MEDCouplingExtrudedMesh(const MEDCouplingUMesh *mesh3D, const MEDCouplingUMesh *mesh2D, int cell2DId) throw(INTERP_KERNEL::Exception)
82 try:_mesh2D(const_cast<MEDCouplingUMesh *>(mesh2D)),_mesh1D(MEDCouplingUMesh::New()),_mesh3D_ids(0),_cell_2D_id(cell2DId)
83 {
84   if(_mesh2D!=0)
85     _mesh2D->incrRef();
86   computeExtrusion(mesh3D);
87   setName(mesh3D->getName());
88 }
89 catch(INTERP_KERNEL::Exception& e)
90   {
91     if(_mesh2D)
92       _mesh2D->decrRef();
93     if(_mesh1D)
94       _mesh1D->decrRef();
95     if(_mesh3D_ids)
96       _mesh3D_ids->decrRef();
97     throw e;
98   }
99
100 MEDCouplingExtrudedMesh::MEDCouplingExtrudedMesh():_mesh2D(0),_mesh1D(0),_mesh3D_ids(0),_cell_2D_id(-1)
101 {
102 }
103
104 MEDCouplingExtrudedMesh::MEDCouplingExtrudedMesh(const MEDCouplingExtrudedMesh& other, bool deepCopy):MEDCouplingMesh(other),_cell_2D_id(other._cell_2D_id)
105 {
106   if(deepCopy)
107     {
108       _mesh2D=other._mesh2D->clone(true);
109       _mesh1D=other._mesh1D->clone(true);
110       _mesh3D_ids=other._mesh3D_ids->deepCpy();
111     }
112   else
113     {
114       _mesh2D=other._mesh2D;
115       if(_mesh2D)
116         _mesh2D->incrRef();
117       _mesh1D=other._mesh1D;
118       if(_mesh1D)
119         _mesh1D->incrRef();
120       _mesh3D_ids=other._mesh3D_ids;
121       if(_mesh3D_ids)
122         _mesh3D_ids->incrRef();
123     }
124 }
125
126 int MEDCouplingExtrudedMesh::getNumberOfCells() const
127 {
128   return _mesh2D->getNumberOfCells()*_mesh1D->getNumberOfCells();
129 }
130
131 int MEDCouplingExtrudedMesh::getNumberOfNodes() const
132 {
133   return _mesh2D->getNumberOfNodes();
134 }
135
136 int MEDCouplingExtrudedMesh::getSpaceDimension() const
137 {
138   return 3;
139 }
140
141 int MEDCouplingExtrudedMesh::getMeshDimension() const
142 {
143   return 3;
144 }
145
146 MEDCouplingMesh *MEDCouplingExtrudedMesh::deepCpy() const
147 {
148   return clone(true);
149 }
150
151 MEDCouplingExtrudedMesh *MEDCouplingExtrudedMesh::clone(bool recDeepCpy) const
152 {
153   return new MEDCouplingExtrudedMesh(*this,recDeepCpy);
154 }
155
156 bool MEDCouplingExtrudedMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
157 {
158   if(!other)
159     throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::isEqualIfNotWhy : input other pointer is null !");
160   const MEDCouplingExtrudedMesh *otherC=dynamic_cast<const MEDCouplingExtrudedMesh *>(other);
161   std::ostringstream oss;
162   if(!otherC)
163     {
164       reason="mesh given in input is not castable in MEDCouplingExtrudedMesh !";
165       return false;
166     }
167   if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
168     return false;
169   if(!_mesh2D->isEqualIfNotWhy(otherC->_mesh2D,prec,reason))
170     {
171       reason.insert(0,"Mesh2D unstructured meshes differ : ");
172       return false;
173     }
174   if(!_mesh1D->isEqualIfNotWhy(otherC->_mesh1D,prec,reason))
175     {
176       reason.insert(0,"Mesh1D unstructured meshes differ : ");
177       return false;
178     }
179   if(!_mesh3D_ids->isEqualIfNotWhy(*otherC->_mesh3D_ids,reason))
180     {
181       reason.insert(0,"Mesh3D ids DataArrayInt instances differ : ");
182       return false;
183     }
184   if(_cell_2D_id!=otherC->_cell_2D_id)
185     {
186       oss << "Cell 2D id of the two extruded mesh differ : this = " << _cell_2D_id << " other = " <<  otherC->_cell_2D_id;
187       reason=oss.str();
188       return false;
189     }
190   return true;
191 }
192
193 bool MEDCouplingExtrudedMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
194 {
195   const MEDCouplingExtrudedMesh *otherC=dynamic_cast<const MEDCouplingExtrudedMesh *>(other);
196   if(!otherC)
197     return false;
198   if(!_mesh2D->isEqualWithoutConsideringStr(otherC->_mesh2D,prec))
199     return false;
200   if(!_mesh1D->isEqualWithoutConsideringStr(otherC->_mesh1D,prec))
201     return false;
202   if(!_mesh3D_ids->isEqualWithoutConsideringStr(*otherC->_mesh3D_ids))
203     return false;
204   if(_cell_2D_id!=otherC->_cell_2D_id)
205     return false;
206   return true;
207 }
208
209 void MEDCouplingExtrudedMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
210                                                    DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
211 {
212   throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::checkDeepEquivalWith : not implemented yet !");
213 }
214
215 void MEDCouplingExtrudedMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
216                                                               DataArrayInt *&cellCor) const throw(INTERP_KERNEL::Exception)
217 {
218   throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::checkDeepEquivalOnSameNodesWith : not implemented yet !");
219 }
220
221 INTERP_KERNEL::NormalizedCellType MEDCouplingExtrudedMesh::getTypeOfCell(int cellId) const
222 {
223   const int *ids=_mesh3D_ids->getConstPointer();
224   int nbOf3DCells=_mesh3D_ids->getNumberOfTuples();
225   const int *where=std::find(ids,ids+nbOf3DCells,cellId);
226   if(where==ids+nbOf3DCells)
227     throw INTERP_KERNEL::Exception("Invalid cellId specified >= getNumberOfCells() !");
228   int nbOfCells2D=_mesh2D->getNumberOfCells();
229   int locId=((int)std::distance(ids,where))%nbOfCells2D;
230   INTERP_KERNEL::NormalizedCellType tmp=_mesh2D->getTypeOfCell(locId);
231   return INTERP_KERNEL::CellModel::GetCellModel(tmp).getExtrudedType();
232 }
233
234 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingExtrudedMesh::getAllGeoTypes() const
235 {
236   const std::set<INTERP_KERNEL::NormalizedCellType>& ret2D=_mesh2D->getAllTypes();
237   std::set<INTERP_KERNEL::NormalizedCellType> ret;
238   for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=ret2D.begin();it!=ret2D.end();it++)
239     ret.insert(INTERP_KERNEL::CellModel::GetCellModel(*it).getExtrudedType());
240   return ret;
241 }
242
243 int MEDCouplingExtrudedMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
244 {
245   int ret=0;
246   int nbOfCells2D=_mesh2D->getNumberOfCells();
247   for(int i=0;i<nbOfCells2D;i++)
248     {
249       INTERP_KERNEL::NormalizedCellType t=_mesh2D->getTypeOfCell(i);
250       if(INTERP_KERNEL::CellModel::GetCellModel(t).getExtrudedType()==type)
251         ret++;
252     }
253   return ret*_mesh1D->getNumberOfCells();
254 }
255
256 void MEDCouplingExtrudedMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
257 {
258   int nbOfCells2D=_mesh2D->getNumberOfCells();
259   int nbOfNodes2D=_mesh2D->getNumberOfNodes();
260   int locId=cellId%nbOfCells2D;
261   int lev=cellId/nbOfCells2D;
262   std::vector<int> tmp,tmp2;
263   _mesh2D->getNodeIdsOfCell(locId,tmp);
264   tmp2=tmp;
265   std::transform(tmp.begin(),tmp.end(),tmp.begin(),std::bind2nd(std::plus<int>(),nbOfNodes2D*lev));
266   std::transform(tmp2.begin(),tmp2.end(),tmp2.begin(),std::bind2nd(std::plus<int>(),nbOfNodes2D*(lev+1)));
267   conn.insert(conn.end(),tmp.begin(),tmp.end());
268   conn.insert(conn.end(),tmp2.begin(),tmp2.end());
269 }
270
271 void MEDCouplingExtrudedMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
272 {
273   int nbOfNodes2D=_mesh2D->getNumberOfNodes();
274   int locId=nodeId%nbOfNodes2D;
275   int lev=nodeId/nbOfNodes2D;
276   std::vector<double> tmp,tmp2;
277   _mesh2D->getCoordinatesOfNode(locId,tmp);
278   tmp2=tmp;
279   int spaceDim=_mesh1D->getSpaceDimension();
280   const double *z=_mesh1D->getCoords()->getConstPointer();
281   std::transform(tmp.begin(),tmp.end(),z+lev*spaceDim,tmp.begin(),std::plus<double>());
282   std::transform(tmp2.begin(),tmp2.end(),z+(lev+1)*spaceDim,tmp2.begin(),std::plus<double>());
283   coo.insert(coo.end(),tmp.begin(),tmp.end());
284   coo.insert(coo.end(),tmp2.begin(),tmp2.end());
285 }
286
287 std::string MEDCouplingExtrudedMesh::simpleRepr() const
288 {
289   std::ostringstream ret;
290   ret << "3D Extruded mesh from a 2D Surf Mesh with name : \"" << getName() << "\"\n";
291   ret << "Description of mesh : \"" << getDescription() << "\"\n";
292   int tmpp1,tmpp2;
293   double tt=getTime(tmpp1,tmpp2);
294   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
295   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
296   ret << "Cell id where 1D mesh has been deduced : " << _cell_2D_id << "\n";
297   ret << "Number of cells : " << getNumberOfCells() << "(" << _mesh2D->getNumberOfCells() << "x" << _mesh1D->getNumberOfCells() << ")\n";
298   ret << "1D Mesh info : _____________________\n\n\n";
299   ret << _mesh1D->simpleRepr();
300   ret << "\n\n\n2D Mesh info : _____________________\n\n\n" << _mesh2D->simpleRepr() << "\n\n\n";
301   return ret.str();
302 }
303
304 std::string MEDCouplingExtrudedMesh::advancedRepr() const
305 {
306   std::ostringstream ret;
307   ret << "3D Extruded mesh from a 2D Surf Mesh with name : \"" << getName() << "\"\n";
308   ret << "Description of mesh : \"" << getDescription() << "\"\n";
309   int tmpp1,tmpp2;
310   double tt=getTime(tmpp1,tmpp2);
311   ret << "Time attached to the mesh (unit) : " << tt << " (" << getTimeUnit() << ")\n";
312   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
313   ret << "Cell id where 1D mesh has been deduced : " << _cell_2D_id << "\n";
314   ret << "Number of cells : " << getNumberOfCells() << "(" << _mesh2D->getNumberOfCells() << "x" << _mesh1D->getNumberOfCells() << ")\n";
315   ret << "1D Mesh info : _____________________\n\n\n";
316   ret << _mesh1D->advancedRepr();
317   ret << "\n\n\n2D Mesh info : _____________________\n\n\n" << _mesh2D->advancedRepr() << "\n\n\n";
318   ret << "3D cell ids per level :\n";
319   return ret.str();
320 }
321
322 void MEDCouplingExtrudedMesh::checkCoherency() const throw (INTERP_KERNEL::Exception)
323 {
324 }
325
326 void MEDCouplingExtrudedMesh::checkCoherency1(double eps) const throw(INTERP_KERNEL::Exception)
327 {
328   checkCoherency();
329 }
330
331 void MEDCouplingExtrudedMesh::checkCoherency2(double eps) const throw(INTERP_KERNEL::Exception)
332 {
333   checkCoherency1(eps);
334 }
335
336 void MEDCouplingExtrudedMesh::getBoundingBox(double *bbox) const
337 {
338   double bbox2D[6];
339   _mesh2D->getBoundingBox(bbox2D);
340   const double *nodes1D=_mesh1D->getCoords()->getConstPointer();
341   int nbOfNodes1D=_mesh1D->getNumberOfNodes();
342   double bbox1DMin[3],bbox1DMax[3],tmp[3];
343   std::fill(bbox1DMin,bbox1DMin+3,std::numeric_limits<double>::max());
344   std::fill(bbox1DMax,bbox1DMax+3,-(std::numeric_limits<double>::max()));
345   for(int i=0;i<nbOfNodes1D;i++)
346     {
347       std::transform(nodes1D+3*i,nodes1D+3*(i+1),bbox1DMin,bbox1DMin,static_cast<const double& (*)(const double&, const double&)>(std::min<double>));
348       std::transform(nodes1D+3*i,nodes1D+3*(i+1),bbox1DMax,bbox1DMax,static_cast<const double& (*)(const double&, const double&)>(std::max<double>));
349     }
350   std::transform(bbox1DMax,bbox1DMax+3,bbox1DMin,tmp,std::minus<double>());
351   int id=(int)std::distance(tmp,std::max_element(tmp,tmp+3));
352   bbox[0]=bbox1DMin[0]; bbox[1]=bbox1DMax[0];
353   bbox[2]=bbox1DMin[1]; bbox[3]=bbox1DMax[1];
354   bbox[4]=bbox1DMin[2]; bbox[5]=bbox1DMax[2];
355   bbox[2*id+1]+=tmp[id];
356 }
357
358 void MEDCouplingExtrudedMesh::updateTime() const
359 {
360   if(_mesh2D)
361     {
362       updateTimeWith(*_mesh2D);
363     }
364   if(_mesh1D)
365     {
366       updateTimeWith(*_mesh1D);
367     }
368 }
369
370 void MEDCouplingExtrudedMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
371 {
372   throw INTERP_KERNEL::Exception("Functionnality of renumbering cells unavailable for ExtrudedMesh");
373 }
374
375 MEDCouplingUMesh *MEDCouplingExtrudedMesh::build3DUnstructuredMesh() const
376 {
377   MEDCouplingUMesh *ret=_mesh2D->buildExtrudedMesh(_mesh1D,0);
378   const int *renum=_mesh3D_ids->getConstPointer();
379   ret->renumberCells(renum,false);
380   ret->setName(getName());
381   return ret;
382 }
383
384 MEDCouplingUMesh *MEDCouplingExtrudedMesh::buildUnstructured() const throw(INTERP_KERNEL::Exception)
385 {
386   return build3DUnstructuredMesh();
387 }
388
389 MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::getMeasureField(bool) const
390 {
391   std::string name="MeasureOfMesh_";
392   name+=getName();
393   MEDCouplingFieldDouble *ret2D=_mesh2D->getMeasureField(true);
394   MEDCouplingFieldDouble *ret1D=_mesh1D->getMeasureField(true);
395   const double *ret2DPtr=ret2D->getArray()->getConstPointer();
396   const double *ret1DPtr=ret1D->getArray()->getConstPointer();
397   int nbOf2DCells=_mesh2D->getNumberOfCells();
398   int nbOf1DCells=_mesh1D->getNumberOfCells();
399   int nbOf3DCells=nbOf2DCells*nbOf1DCells;
400   const int *renum=_mesh3D_ids->getConstPointer();
401   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
402   ret->setMesh(this);
403   DataArrayDouble *da=DataArrayDouble::New();
404   da->alloc(nbOf3DCells,1);
405   double *retPtr=da->getPointer();
406   for(int i=0;i<nbOf1DCells;i++)
407     for(int j=0;j<nbOf2DCells;j++)
408       retPtr[renum[i*nbOf2DCells+j]]=ret2DPtr[j]*ret1DPtr[i];
409   ret->setArray(da);
410   da->decrRef();
411   ret->setName(name.c_str());
412   ret2D->decrRef();
413   ret1D->decrRef();
414   return ret;
415 }
416
417 MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::getMeasureFieldOnNode(bool isAbs) const
418 {
419   //not implemented yet
420   return 0;
421 }
422
423 MEDCouplingFieldDouble *MEDCouplingExtrudedMesh::buildOrthogonalField() const
424 {
425   throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::buildOrthogonalField : This method has no sense for MEDCouplingExtrudedMesh that is 3D !");
426 }
427
428 int MEDCouplingExtrudedMesh::getCellContainingPoint(const double *pos, double eps) const
429 {
430   throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::getCellContainingPoint : not implemented yet !");
431 }
432
433 MEDCouplingExtrudedMesh::~MEDCouplingExtrudedMesh()
434 {
435   if(_mesh2D)
436     _mesh2D->decrRef();
437   if(_mesh1D)
438     _mesh1D->decrRef();
439   if(_mesh3D_ids)
440     _mesh3D_ids->decrRef();
441 }
442
443 void MEDCouplingExtrudedMesh::computeExtrusion(const MEDCouplingUMesh *mesh3D) throw(INTERP_KERNEL::Exception)
444 {
445   const char errMsg1[]="2D mesh is empty unable to compute extrusion !";
446   const char errMsg2[]="Coords between 2D and 3D meshes are not the same ! Try MEDCouplingPointSet::tryToShareSameCoords method";
447   const char errMsg3[]="No chance to find extrusion pattern in mesh3D,mesh2D couple because nbCells3D%nbCells2D!=0 !";
448   if(_mesh2D==0 || mesh3D==0)
449     throw INTERP_KERNEL::Exception(errMsg1);
450   if(_mesh2D->getCoords()!=mesh3D->getCoords())
451     throw INTERP_KERNEL::Exception(errMsg2);
452   if(mesh3D->getNumberOfCells()%_mesh2D->getNumberOfCells()!=0)
453     throw INTERP_KERNEL::Exception(errMsg3);
454   if(!_mesh3D_ids)
455     _mesh3D_ids=DataArrayInt::New();
456   if(!_mesh1D)
457     _mesh1D=MEDCouplingUMesh::New();
458   computeExtrusionAlg(mesh3D);
459 }
460
461 void MEDCouplingExtrudedMesh::build1DExtrusion(int idIn3DDesc, int newId, int nbOf1DLev, MEDCouplingUMesh *subMesh,
462                                                const int *desc3D, const int *descIndx3D,
463                                                const int *revDesc3D, const int *revDescIndx3D,
464                                                bool computeMesh1D) throw(INTERP_KERNEL::Exception)
465 {
466   int nbOf2DCells=_mesh2D->getNumberOfCells();
467   int start=revDescIndx3D[idIn3DDesc];
468   int end=revDescIndx3D[idIn3DDesc+1];
469   if(end-start!=1)
470     {
471       std::ostringstream ost; ost << "Invalid bases 2D mesh specified : 2D cell # " <<  idIn3DDesc;
472       ost << " shared by more than 1 3D cell !!!";
473       throw INTERP_KERNEL::Exception(ost.str().c_str());
474     }
475   int current3DCell=revDesc3D[start];
476   int current2DCell=idIn3DDesc;
477   int *mesh3DIDs=_mesh3D_ids->getPointer();
478   mesh3DIDs[newId]=current3DCell;
479   const int *conn2D=subMesh->getNodalConnectivity()->getConstPointer();
480   const int *conn2DIndx=subMesh->getNodalConnectivityIndex()->getConstPointer();
481   for(int i=1;i<nbOf1DLev;i++)
482     {
483       std::vector<int> conn(conn2D+conn2DIndx[current2DCell]+1,conn2D+conn2DIndx[current2DCell+1]);
484       std::sort(conn.begin(),conn.end());
485       if(computeMesh1D)
486         computeBaryCenterOfFace(conn,i-1);
487       current2DCell=findOppositeFaceOf(current2DCell,current3DCell,conn,
488                                        desc3D,descIndx3D,conn2D,conn2DIndx);
489       start=revDescIndx3D[current2DCell];
490       end=revDescIndx3D[current2DCell+1];
491       if(end-start!=2)
492         {
493           std::ostringstream ost; ost << "Expecting to have 2 3D cells attached to 2D cell " << current2DCell << "!";
494           ost << " : Impossible or call tryToShareSameCoords method !";
495           throw INTERP_KERNEL::Exception(ost.str().c_str());
496         }
497       if(revDesc3D[start]!=current3DCell)
498         current3DCell=revDesc3D[start];
499       else
500         current3DCell=revDesc3D[start+1];
501       mesh3DIDs[i*nbOf2DCells+newId]=current3DCell;
502     }
503   if(computeMesh1D)
504     {
505       std::vector<int> conn(conn2D+conn2DIndx[current2DCell]+1,conn2D+conn2DIndx[current2DCell+1]);
506       std::sort(conn.begin(),conn.end());
507       computeBaryCenterOfFace(conn,nbOf1DLev-1);
508       current2DCell=findOppositeFaceOf(current2DCell,current3DCell,conn,
509                                        desc3D,descIndx3D,conn2D,conn2DIndx);
510       conn.clear();
511       conn.insert(conn.end(),conn2D+conn2DIndx[current2DCell]+1,conn2D+conn2DIndx[current2DCell+1]);
512       std::sort(conn.begin(),conn.end());
513       computeBaryCenterOfFace(conn,nbOf1DLev);
514     }
515 }
516
517 int MEDCouplingExtrudedMesh::findOppositeFaceOf(int current2DCell, int current3DCell, const std::vector<int>& connSorted,
518                                                 const int *desc3D, const int *descIndx3D,
519                                                 const int *conn2D, const int *conn2DIndx) throw(INTERP_KERNEL::Exception)
520 {
521   int start=descIndx3D[current3DCell];
522   int end=descIndx3D[current3DCell+1];
523   bool found=false;
524   for(const int *candidate2D=desc3D+start;candidate2D!=desc3D+end && !found;candidate2D++)
525     {
526       if(*candidate2D!=current2DCell)
527         {
528           std::vector<int> conn2(conn2D+conn2DIndx[*candidate2D]+1,conn2D+conn2DIndx[*candidate2D+1]);
529           std::sort(conn2.begin(),conn2.end());
530           std::list<int> intersect;
531           std::set_intersection(connSorted.begin(),connSorted.end(),conn2.begin(),conn2.end(),
532                                 std::insert_iterator< std::list<int> >(intersect,intersect.begin()));
533           if(intersect.empty())
534             return *candidate2D;
535         }
536     }
537   std::ostringstream ost; ost << "Impossible to find an opposite 2D face of face # " <<  current2DCell;
538   ost << " in 3D cell # " << current3DCell << " : Impossible or call tryToShareSameCoords method !";
539   throw INTERP_KERNEL::Exception(ost.str().c_str());
540 }
541
542 void MEDCouplingExtrudedMesh::computeBaryCenterOfFace(const std::vector<int>& nodalConnec, int lev1DId)
543 {
544   double *zoneToUpdate=_mesh1D->getCoords()->getPointer()+lev1DId*3;
545   std::fill(zoneToUpdate,zoneToUpdate+3,0.);
546   const double *coords=_mesh2D->getCoords()->getConstPointer();
547   for(std::vector<int>::const_iterator iter=nodalConnec.begin();iter!=nodalConnec.end();iter++)
548     std::transform(zoneToUpdate,zoneToUpdate+3,coords+3*(*iter),zoneToUpdate,std::plus<double>());
549   std::transform(zoneToUpdate,zoneToUpdate+3,zoneToUpdate,std::bind2nd(std::multiplies<double>(),(double)(1./(int)nodalConnec.size())));
550 }
551
552 int MEDCouplingExtrudedMesh::FindCorrespCellByNodalConn(const std::vector<int>& nodalConnec, const int *revNodalPtr, const int *revNodalIndxPtr) throw(INTERP_KERNEL::Exception)
553 {
554   std::vector<int>::const_iterator iter=nodalConnec.begin();
555   std::set<int> s1(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1]);
556   iter++;
557   for(;iter!=nodalConnec.end();iter++)
558     {
559       std::set<int> s2(revNodalPtr+revNodalIndxPtr[*iter],revNodalPtr+revNodalIndxPtr[*iter+1]);
560       std::set<int> s3;
561       std::set_intersection(s1.begin(),s1.end(),s2.begin(),s2.end(),std::insert_iterator< std::set<int> >(s3,s3.end()));
562       s1=s3;
563     }
564   if(s1.size()==1)
565     return *(s1.begin());
566   std::ostringstream ostr;
567   ostr << "Cell with nodal connec : ";
568   std::copy(nodalConnec.begin(),nodalConnec.end(),std::ostream_iterator<int>(ostr," "));
569   ostr << " is not part of mesh";
570   throw INTERP_KERNEL::Exception(ostr.str().c_str());
571 }
572
573 /*!
574  * This method is callable on 1Dmeshes (meshDim==1 && spaceDim==3) returned by MEDCouplingExtrudedMesh::getMesh1D typically.
575  * 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.
576  * 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
577  * are created ('m1r' and 'm2r') that can be used for interpolation.
578  * @param m1 input mesh with meshDim==1 and spaceDim==3
579  * @param m2 input mesh with meshDim==1 and spaceDim==3
580  * @param eps tolerance acceptable to determine compatibility
581  * @param m1r output mesh with ref count equal to 1 with meshDim==1 and spaceDim==1
582  * @param m2r output mesh with ref count equal to 1 with meshDim==1 and spaceDim==1
583  * @param v is the output normalized vector of the common direction of 'm1' and 'm2'  
584  * @throw in case that m1 and m2 are not compatible each other.
585  */
586 void MEDCouplingExtrudedMesh::Project1DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps,
587                                               MEDCouplingUMesh *&m1r, MEDCouplingUMesh *&m2r, double *v) throw(INTERP_KERNEL::Exception)
588 {
589   if(m1->getSpaceDimension()!=3 || m1->getSpaceDimension()!=3)
590     throw INTERP_KERNEL::Exception("Input meshes are expected to have a spaceDim==3 for Projec1D !");
591   m1r=m1->clone(true);
592   m2r=m2->clone(true);
593   m1r->changeSpaceDimension(1);
594   m2r->changeSpaceDimension(1);
595   std::vector<int> c;
596   std::vector<double> ref,ref2;
597   m1->getNodeIdsOfCell(0,c);
598   m1->getCoordinatesOfNode(c[0],ref);
599   m1->getCoordinatesOfNode(c[1],ref2);
600   std::transform(ref2.begin(),ref2.end(),ref.begin(),v,std::minus<double>());
601   double n=INTERP_KERNEL::norm<3>(v);
602   std::transform(v,v+3,v,std::bind2nd(std::multiplies<double>(),1/n));
603   m1->project1D(&ref[0],v,eps,m1r->getCoords()->getPointer());
604   m2->project1D(&ref[0],v,eps,m2r->getCoords()->getPointer());
605   
606 }
607
608 void MEDCouplingExtrudedMesh::rotate(const double *center, const double *vector, double angle)
609 {
610   _mesh2D->rotate(center,vector,angle);
611   _mesh1D->rotate(center,vector,angle);
612 }
613
614 void MEDCouplingExtrudedMesh::translate(const double *vector)
615 {
616   _mesh2D->translate(vector);
617   _mesh1D->translate(vector);
618 }
619
620 void MEDCouplingExtrudedMesh::scale(const double *point, double factor)
621 {
622   _mesh2D->scale(point,factor);
623   _mesh1D->scale(point,factor);
624 }
625
626 std::vector<int> MEDCouplingExtrudedMesh::getDistributionOfTypes() const throw(INTERP_KERNEL::Exception)
627 {
628   throw INTERP_KERNEL::Exception("Not implemented yet !");
629 }
630
631 DataArrayInt *MEDCouplingExtrudedMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
632 {
633   throw INTERP_KERNEL::Exception("Not implemented yet !");
634 }
635
636 void MEDCouplingExtrudedMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
637 {
638   throw INTERP_KERNEL::Exception("Not implemented yet !");
639 }
640
641 MEDCouplingMesh *MEDCouplingExtrudedMesh::buildPart(const int *start, const int *end) const
642 {
643   // not implemented yet !
644   return 0;
645 }
646
647 MEDCouplingMesh *MEDCouplingExtrudedMesh::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
648 {
649   // not implemented yet !
650   return 0;
651 }
652
653 DataArrayInt *MEDCouplingExtrudedMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception)
654 {
655   throw INTERP_KERNEL::Exception("MEDCouplingExtrudedMesh::simplexize : unavailable for such a type of mesh : Extruded !");
656 }
657
658 MEDCouplingMesh *MEDCouplingExtrudedMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
659 {
660   // not implemented yet !
661   return 0;
662 }
663
664 DataArrayDouble *MEDCouplingExtrudedMesh::getCoordinatesAndOwner() const
665 {
666   DataArrayDouble *arr2D=_mesh2D->getCoords();
667   DataArrayDouble *arr1D=_mesh1D->getCoords();
668   DataArrayDouble *ret=DataArrayDouble::New();
669   ret->alloc(getNumberOfNodes(),3);
670   int nbOf1DLev=_mesh1D->getNumberOfNodes();
671   int nbOf2DNodes=_mesh2D->getNumberOfNodes();
672   const double *ptSrc=arr2D->getConstPointer();
673   double *pt=ret->getPointer();
674   std::copy(ptSrc,ptSrc+3*nbOf2DNodes,pt);
675   for(int i=1;i<nbOf1DLev;i++)
676     {
677       std::copy(ptSrc,ptSrc+3*nbOf2DNodes,pt+3*i*nbOf2DNodes);
678       double vec[3];
679       std::copy(arr1D->getConstPointer()+3*i,arr1D->getConstPointer()+3*(i+1),vec);
680       std::transform(arr1D->getConstPointer()+3*(i-1),arr1D->getConstPointer()+3*i,vec,vec,std::minus<double>());
681       for(int j=0;j<nbOf2DNodes;j++)
682         std::transform(vec,vec+3,pt+3*(i*nbOf2DNodes+j),pt+3*(i*nbOf2DNodes+j),std::plus<double>());
683     }
684   return ret;
685 }
686
687 DataArrayDouble *MEDCouplingExtrudedMesh::getBarycenterAndOwner() const
688 {
689   //not yet implemented
690   return 0;
691 }
692
693 void MEDCouplingExtrudedMesh::computeExtrusionAlg(const MEDCouplingUMesh *mesh3D) throw(INTERP_KERNEL::Exception)
694 {
695   _mesh3D_ids->alloc(mesh3D->getNumberOfCells(),1);
696   int nbOf1DLev=mesh3D->getNumberOfCells()/_mesh2D->getNumberOfCells();
697   _mesh1D->setMeshDimension(1);
698   _mesh1D->allocateCells(nbOf1DLev);
699   int tmpConn[2];
700   for(int i=0;i<nbOf1DLev;i++)
701     {
702       tmpConn[0]=i;
703       tmpConn[1]=i+1;
704       _mesh1D->insertNextCell(INTERP_KERNEL::NORM_SEG2,2,tmpConn);
705     }
706   _mesh1D->finishInsertingCells();
707   DataArrayDouble *myCoords=DataArrayDouble::New();
708   myCoords->alloc(nbOf1DLev+1,3);
709   _mesh1D->setCoords(myCoords);
710   myCoords->decrRef();
711   DataArrayInt *desc,*descIndx,*revDesc,*revDescIndx;
712   desc=DataArrayInt::New(); descIndx=DataArrayInt::New(); revDesc=DataArrayInt::New(); revDescIndx=DataArrayInt::New();
713   MEDCouplingUMesh *subMesh=mesh3D->buildDescendingConnectivity(desc,descIndx,revDesc,revDescIndx);
714   DataArrayInt *revNodal2D,*revNodalIndx2D;
715   revNodal2D=DataArrayInt::New(); revNodalIndx2D=DataArrayInt::New();
716   subMesh->getReverseNodalConnectivity(revNodal2D,revNodalIndx2D);
717   const int *nodal2D=_mesh2D->getNodalConnectivity()->getConstPointer();
718   const int *nodal2DIndx=_mesh2D->getNodalConnectivityIndex()->getConstPointer();
719   const int *revNodal2DPtr=revNodal2D->getConstPointer();
720   const int *revNodalIndx2DPtr=revNodalIndx2D->getConstPointer();
721   const int *descP=desc->getConstPointer();
722   const int *descIndxP=descIndx->getConstPointer();
723   const int *revDescP=revDesc->getConstPointer();
724   const int *revDescIndxP=revDescIndx->getConstPointer();
725   //
726   int nbOf2DCells=_mesh2D->getNumberOfCells();
727   for(int i=0;i<nbOf2DCells;i++)
728     {
729       int idInSubMesh;
730        std::vector<int> nodalConnec(nodal2D+nodal2DIndx[i]+1,nodal2D+nodal2DIndx[i+1]);
731        try
732         {
733           idInSubMesh=FindCorrespCellByNodalConn(nodalConnec,revNodal2DPtr,revNodalIndx2DPtr);
734         }
735        catch(INTERP_KERNEL::Exception& e)
736          {
737            std::ostringstream ostr; ostr << "mesh2D cell # " << i << " is not part of any cell of 3D mesh !\n";
738            ostr << e.what();
739            throw INTERP_KERNEL::Exception(ostr.str().c_str());
740          }
741        build1DExtrusion(idInSubMesh,i,nbOf1DLev,subMesh,descP,descIndxP,revDescP,revDescIndxP,i==_cell_2D_id);
742     }
743   //
744   revNodal2D->decrRef();
745   revNodalIndx2D->decrRef();
746   subMesh->decrRef();
747   desc->decrRef();
748   descIndx->decrRef();
749   revDesc->decrRef();
750   revDescIndx->decrRef();
751 }
752
753 void MEDCouplingExtrudedMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
754 {
755   std::vector<int> tinyInfo1;
756   std::vector<std::string> ls1;
757   std::vector<double> ls3;
758   _mesh2D->getTinySerializationInformation(ls3,tinyInfo1,ls1);
759   std::vector<int> tinyInfo2;
760   std::vector<std::string> ls2;
761   std::vector<double> ls4;
762   _mesh1D->getTinySerializationInformation(ls4,tinyInfo2,ls2);
763   tinyInfo.clear(); littleStrings.clear();
764   tinyInfo.insert(tinyInfo.end(),tinyInfo1.begin(),tinyInfo1.end());
765   littleStrings.insert(littleStrings.end(),ls1.begin(),ls1.end());
766   tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
767   littleStrings.insert(littleStrings.end(),ls2.begin(),ls2.end());
768   tinyInfo.push_back(_cell_2D_id);
769   tinyInfo.push_back((int)tinyInfo1.size());
770   tinyInfo.push_back(_mesh3D_ids->getNbOfElems());
771   littleStrings.push_back(getName());
772   littleStrings.push_back(getDescription());
773 }
774
775 void MEDCouplingExtrudedMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
776 {
777   std::size_t sz=tinyInfo.size();
778   int sz1=tinyInfo[sz-2];
779   std::vector<int> ti1(tinyInfo.begin(),tinyInfo.begin()+sz1);
780   std::vector<int> ti2(tinyInfo.begin()+sz1,tinyInfo.end()-3);
781   MEDCouplingUMesh *um=MEDCouplingUMesh::New();
782   DataArrayInt *a1tmp=DataArrayInt::New();
783   DataArrayDouble *a2tmp=DataArrayDouble::New();
784   int la1=0,la2=0;
785   std::vector<std::string> ls1,ls2;
786   um->resizeForUnserialization(ti1,a1tmp,a2tmp,ls1);
787   la1+=a1tmp->getNbOfElems(); la2+=a2tmp->getNbOfElems();
788   a1tmp->decrRef(); a2tmp->decrRef();
789   a1tmp=DataArrayInt::New(); a2tmp=DataArrayDouble::New();
790   um->resizeForUnserialization(ti2,a1tmp,a2tmp,ls2);
791   la1+=a1tmp->getNbOfElems(); la2+=a2tmp->getNbOfElems();
792   a1tmp->decrRef(); a2tmp->decrRef();
793   um->decrRef();
794   //
795   a1->alloc(la1+tinyInfo[sz-1],1);
796   a2->alloc(la2,1);
797   littleStrings.resize(ls1.size()+ls2.size()+2);
798 }
799
800 void MEDCouplingExtrudedMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
801 {
802   a1=DataArrayInt::New(); a2=DataArrayDouble::New();
803   DataArrayInt *a1_1=0,*a1_2=0;
804   DataArrayDouble *a2_1=0,*a2_2=0;
805   _mesh2D->serialize(a1_1,a2_1);
806   _mesh1D->serialize(a1_2,a2_2);
807   a1->alloc(a1_1->getNbOfElems()+a1_2->getNbOfElems()+_mesh3D_ids->getNbOfElems(),1);
808   int *ptri=a1->getPointer();
809   ptri=std::copy(a1_1->getConstPointer(),a1_1->getConstPointer()+a1_1->getNbOfElems(),ptri);
810   a1_1->decrRef();
811   ptri=std::copy(a1_2->getConstPointer(),a1_2->getConstPointer()+a1_2->getNbOfElems(),ptri);
812   a1_2->decrRef();
813   std::copy(_mesh3D_ids->getConstPointer(),_mesh3D_ids->getConstPointer()+_mesh3D_ids->getNbOfElems(),ptri);
814   a2->alloc(a2_1->getNbOfElems()+a2_2->getNbOfElems(),1);
815   double *ptrd=a2->getPointer();
816   ptrd=std::copy(a2_1->getConstPointer(),a2_1->getConstPointer()+a2_1->getNbOfElems(),ptrd);
817   a2_1->decrRef();
818   std::copy(a2_2->getConstPointer(),a2_2->getConstPointer()+a2_2->getNbOfElems(),ptrd);
819   a2_2->decrRef();
820 }
821
822 void MEDCouplingExtrudedMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2, const std::vector<std::string>& littleStrings)
823 {
824   setName(littleStrings[littleStrings.size()-2].c_str());
825   setDescription(littleStrings.back().c_str());
826   std::size_t sz=tinyInfo.size();
827   int sz1=tinyInfo[sz-2];
828   _cell_2D_id=tinyInfo[sz-3];
829   std::vector<int> ti1(tinyInfo.begin(),tinyInfo.begin()+sz1);
830   std::vector<int> ti2(tinyInfo.begin()+sz1,tinyInfo.end()-3);
831   DataArrayInt *a1tmp=DataArrayInt::New();
832   DataArrayDouble *a2tmp=DataArrayDouble::New();
833   const int *a1Ptr=a1->getConstPointer();
834   const double *a2Ptr=a2->getConstPointer();
835   _mesh2D=MEDCouplingUMesh::New();
836   std::vector<std::string> ls1,ls2;
837   _mesh2D->resizeForUnserialization(ti1,a1tmp,a2tmp,ls1);
838   std::copy(a2Ptr,a2Ptr+a2tmp->getNbOfElems(),a2tmp->getPointer());
839   std::copy(a1Ptr,a1Ptr+a1tmp->getNbOfElems(),a1tmp->getPointer());
840   a2Ptr+=a2tmp->getNbOfElems();
841   a1Ptr+=a1tmp->getNbOfElems();
842   ls2.insert(ls2.end(),littleStrings.begin(),littleStrings.begin()+ls1.size());
843   std::vector<double> d1(1);
844   _mesh2D->unserialization(d1,ti1,a1tmp,a2tmp,ls2);
845   a1tmp->decrRef(); a2tmp->decrRef();
846   //
847   ls2.clear();
848   ls2.insert(ls2.end(),littleStrings.begin()+ls1.size(),littleStrings.end()-2);
849   _mesh1D=MEDCouplingUMesh::New();
850   a1tmp=DataArrayInt::New(); a2tmp=DataArrayDouble::New();
851   _mesh1D->resizeForUnserialization(ti2,a1tmp,a2tmp,ls1);
852   std::copy(a2Ptr,a2Ptr+a2tmp->getNbOfElems(),a2tmp->getPointer());
853   std::copy(a1Ptr,a1Ptr+a1tmp->getNbOfElems(),a1tmp->getPointer());
854   a1Ptr+=a1tmp->getNbOfElems();
855   _mesh1D->unserialization(d1,ti2,a1tmp,a2tmp,ls2);
856   a1tmp->decrRef(); a2tmp->decrRef();
857   //
858   _mesh3D_ids=DataArrayInt::New();
859   int szIds=(int)std::distance(a1Ptr,a1->getConstPointer()+a1->getNbOfElems());
860   _mesh3D_ids->alloc(szIds,1);
861   std::copy(a1Ptr,a1Ptr+szIds,_mesh3D_ids->getPointer());
862 }
863
864 void MEDCouplingExtrudedMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception)
865 {
866   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> m=buildUnstructured();
867   m->writeVTKLL(ofs,cellData,pointData);
868 }
869
870 std::string MEDCouplingExtrudedMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception)
871 {
872   return _mesh2D->getVTKDataSetType();
873 }