1 // Copyright (C) 2007-2012 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 "MEDCouplingCurveLinearMesh.hxx"
22 #include "MEDCouplingPointSet.hxx"
23 #include "MEDCouplingMemArray.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
26 #include "VolSurfUser.txx"
27 #include "PointLocatorAlgos.txx"
34 using namespace ParaMEDMEM;
36 MEDCouplingCurveLinearMesh::MEDCouplingCurveLinearMesh():_coords(0),_structure(0)
40 MEDCouplingCurveLinearMesh::MEDCouplingCurveLinearMesh(const MEDCouplingCurveLinearMesh& other, bool deepCopy):MEDCouplingStructuredMesh(other,deepCopy),_structure(other._structure)
44 if((const DataArrayDouble *)other._coords)
45 _coords=other._coords->deepCpy();
48 _coords=other._coords;
51 MEDCouplingCurveLinearMesh::~MEDCouplingCurveLinearMesh()
55 MEDCouplingCurveLinearMesh *MEDCouplingCurveLinearMesh::New()
57 return new MEDCouplingCurveLinearMesh;
60 MEDCouplingCurveLinearMesh *MEDCouplingCurveLinearMesh::New(const char *meshName)
62 MEDCouplingCurveLinearMesh *ret=new MEDCouplingCurveLinearMesh;
63 ret->setName(meshName);
67 MEDCouplingMesh *MEDCouplingCurveLinearMesh::deepCpy() const
72 MEDCouplingCurveLinearMesh *MEDCouplingCurveLinearMesh::clone(bool recDeepCpy) const
74 return new MEDCouplingCurveLinearMesh(*this,recDeepCpy);
77 void MEDCouplingCurveLinearMesh::updateTime() const
79 if((const DataArrayDouble *)_coords)
80 updateTimeWith(*_coords);
83 std::size_t MEDCouplingCurveLinearMesh::getHeapMemorySize() const
86 ret+=_structure.capacity()*sizeof(int);
87 if((const DataArrayDouble *)_coords)
88 ret+=_coords->getHeapMemorySize();
89 return MEDCouplingStructuredMesh::getHeapMemorySize()+ret;
93 * This method copyies all tiny strings from other (name and components name).
94 * @throw if other and this have not same mesh type.
96 void MEDCouplingCurveLinearMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
98 const MEDCouplingCurveLinearMesh *otherC=dynamic_cast<const MEDCouplingCurveLinearMesh *>(other);
100 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::copyTinyStringsFrom : meshes have not same type !");
101 MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
102 if((DataArrayDouble *)_coords && (const DataArrayDouble *)otherC->_coords)
103 _coords->copyStringInfoFrom(*otherC->_coords);
106 bool MEDCouplingCurveLinearMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
109 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::isEqualIfNotWhy : input other pointer is null !");
110 const MEDCouplingCurveLinearMesh *otherC=dynamic_cast<const MEDCouplingCurveLinearMesh *>(other);
113 reason="mesh given in input is not castable in MEDCouplingCurveLinearMesh !";
116 if(!MEDCouplingStructuredMesh::isEqualIfNotWhy(other,prec,reason))
118 std::ostringstream oss; oss.precision(15);
119 if(((const DataArrayDouble *)_coords && ((const DataArrayDouble *)otherC->_coords)==0) || (((const DataArrayDouble *)_coords)==0 && (const DataArrayDouble *)otherC->_coords))
121 oss << "Only one CurveLinearMesh between the two this and other has its coordinates defined !";
125 if((const DataArrayDouble *)_coords)
127 if(!_coords->isEqualIfNotWhy(*(otherC->_coords),prec,reason))
129 oss << "Coordinates DataArrayDouble of differ :";
130 reason.insert(0,oss.str());
133 if(_structure!=otherC->_structure)
134 { reason="CurveLinearMesh structures differ !"; return false; }
139 bool MEDCouplingCurveLinearMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
141 const MEDCouplingCurveLinearMesh *otherC=dynamic_cast<const MEDCouplingCurveLinearMesh *>(other);
144 if(((const DataArrayDouble *)_coords && ((const DataArrayDouble *)otherC->_coords)==0) || (((const DataArrayDouble *)_coords)==0 && (const DataArrayDouble *)otherC->_coords))
146 if((const DataArrayDouble *)_coords)
148 if(!_coords->isEqualWithoutConsideringStr(*(otherC->_coords),prec))
150 if(_structure!=otherC->_structure)
156 void MEDCouplingCurveLinearMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
157 DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
159 if(!isEqualWithoutConsideringStr(other,prec))
160 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkDeepEquivalWith : Meshes are not the same !");
164 * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingCurveLinearMesh instance too).
165 * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingCurveLinearMesh, 'this' and 'other' are the same !
167 void MEDCouplingCurveLinearMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
168 DataArrayInt *&cellCor) const throw(INTERP_KERNEL::Exception)
170 const MEDCouplingCurveLinearMesh *otherC=dynamic_cast<const MEDCouplingCurveLinearMesh *>(other);
172 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkDeepEquivalOnSameNodesWith : other is NOT a cartesian mesh ! Impossible to check equivalence !");
175 void MEDCouplingCurveLinearMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
177 std::size_t sz=_structure.size(),i=0,nbOfNodes=1;
179 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkCoherency : structure should have a lgth of size 1 at least !");
180 for(std::vector<int>::const_iterator it=_structure.begin();it!=_structure.end();it++,i++)
183 { std::ostringstream oss; oss << "MEDCouplingCurveLinearMesh::checkCoherency : At pos #" << i << " of structure value is " << *it << "should be >= 1 !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
186 if(!((const DataArrayDouble *)_coords))
187 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkCoherency : the array is not set !");
188 if(!_coords->isAllocated())
189 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkCoherency : the array is not allocated !");
190 if(_coords->getNumberOfComponents()<1)
191 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkCoherency : the array should have >= 1 components !");
192 if(_coords->getNumberOfTuples()!=(int)nbOfNodes)
194 std::ostringstream oss; oss << "MEDCouplingCurveLinearMesh::checkCoherency : structure said that number of nodes should be equal to " << nbOfNodes << " but number of tuples in array is equal to " << _coords->getNumberOfTuples() << " !";
195 throw INTERP_KERNEL::Exception(oss.str().c_str());
199 void MEDCouplingCurveLinearMesh::checkCoherency1(double eps) const throw(INTERP_KERNEL::Exception)
204 void MEDCouplingCurveLinearMesh::checkCoherency2(double eps) const throw(INTERP_KERNEL::Exception)
206 checkCoherency1(eps);
209 int MEDCouplingCurveLinearMesh::getNumberOfCells() const
212 std::size_t nbOfCells=1,i=0;
213 for(std::vector<int>::const_iterator it=_structure.begin();it!=_structure.end();it++,i++)
215 return (int)nbOfCells;
218 int MEDCouplingCurveLinearMesh::getNumberOfNodes() const
221 std::size_t nbOfNodes=1;
222 for(std::vector<int>::const_iterator it=_structure.begin();it!=_structure.end();it++)
224 return (int)nbOfNodes;
227 void MEDCouplingCurveLinearMesh::getSplitCellValues(int *res) const
229 int meshDim=getMeshDimension();
230 for(int l=0;l<meshDim;l++)
233 for(int p=0;p<meshDim-l-1;p++)
234 val*=_structure[p]-1;
235 res[meshDim-l-1]=val;
239 void MEDCouplingCurveLinearMesh::getSplitNodeValues(int *res) const
241 int meshDim=getMeshDimension();
242 for(int l=0;l<meshDim;l++)
245 for(int p=0;p<meshDim-l-1;p++)
247 res[meshDim-l-1]=val;
251 void MEDCouplingCurveLinearMesh::getNodeGridStructure(int *res) const
253 std::copy(_structure.begin(),_structure.end(),res);
256 int MEDCouplingCurveLinearMesh::getSpaceDimension() const
258 if(!((const DataArrayDouble *)_coords))
259 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getSpaceDimension : no array set ! impossible to deduce a space dimension !");
260 return _coords->getNumberOfComponents();
263 int MEDCouplingCurveLinearMesh::getMeshDimension() const
265 return (int)_structure.size();
268 void MEDCouplingCurveLinearMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
270 if(!((const DataArrayDouble *)_coords))
271 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getCoordinatesOfNode : Coordinates not set !");
272 int nbOfCompo=_coords->getNumberOfComponents();
273 if(nodeId>=0 && nodeId<_coords->getNumberOfTuples())
274 coo.insert(coo.end(),_coords->begin()+nodeId*nbOfCompo,_coords->begin()+(nodeId+1)*nbOfCompo);
276 { std::ostringstream oss; oss << "MEDCouplingCurveLinearMesh::getCoordinatesOfNode : nodeId has to be in [0," << _coords->getNumberOfTuples() << ") !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
279 std::string MEDCouplingCurveLinearMesh::simpleRepr() const
281 std::ostringstream ret;
282 ret << "Curve linear mesh with name : \"" << getName() << "\"\n";
283 ret << "Description of mesh : \"" << getDescription() << "\"\n";
285 double tt=getTime(tmpp1,tmpp2);
286 ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
287 ret << "Iteration : " << tmpp1 << " Order : " << tmpp2 << "\n";
288 ret << "The nodal stucture of curve linear mesh is : [";
289 std::copy(_structure.begin(),_structure.end(),std::ostream_iterator<int>(ret,",")); ret << "]\n";
290 ret << "The coords array is this : ";
291 if((const DataArrayDouble *)_coords)
292 _coords->reprZipWithoutNameStream(ret);
294 ret << "no array specified !";
298 std::string MEDCouplingCurveLinearMesh::advancedRepr() const
303 DataArrayDouble *MEDCouplingCurveLinearMesh::getCoords() throw(INTERP_KERNEL::Exception)
308 const DataArrayDouble *MEDCouplingCurveLinearMesh::getCoords() const throw(INTERP_KERNEL::Exception)
313 void MEDCouplingCurveLinearMesh::setCoords(const DataArrayDouble *coords) throw(INTERP_KERNEL::Exception)
315 if(coords!=(const DataArrayDouble *)_coords)
317 _coords=const_cast<DataArrayDouble *>(coords);
324 void MEDCouplingCurveLinearMesh::setNodeGridStructure(const int *gridStructBg, const int *gridStructEnd) throw(INTERP_KERNEL::Exception)
326 _structure.resize(0);
327 _structure.insert(_structure.end(),gridStructBg,gridStructEnd);
330 std::vector<int> MEDCouplingCurveLinearMesh::getNodeGridStructure() const throw(INTERP_KERNEL::Exception)
335 void MEDCouplingCurveLinearMesh::getBoundingBox(double *bbox) const
337 if(!((const DataArrayDouble *)_coords))
338 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getBoundingBox : Coordinates not set !");
339 _coords->getMinMaxPerComponent(bbox);
342 MEDCouplingFieldDouble *MEDCouplingCurveLinearMesh::getMeasureField(bool isAbs) const
345 int meshDim=getMeshDimension();
346 std::string name="MeasureOfMesh_"; name+=getName();
347 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> field=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
348 field->setName(name.c_str()); field->setMesh(const_cast<MEDCouplingCurveLinearMesh *>(this)); field->synchronizeTimeWithMesh();
352 { getMeasureFieldMeshDim3(isAbs,field); return field.retn(); }
354 { getMeasureFieldMeshDim2(isAbs,field); return field.retn(); }
356 { getMeasureFieldMeshDim1(isAbs,field); return field.retn(); }
358 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getMeasureField : mesh dimension must be in [1,2,3] !");
363 * \param [in,out] f field feeded with good values.
364 * \sa MEDCouplingCurveLinearMesh::getMeasureField
366 void MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim1(bool isAbs, MEDCouplingFieldDouble *field) const throw(INTERP_KERNEL::Exception)
368 int nbnodes=getNumberOfNodes();
369 int spaceDim=getSpaceDimension();
370 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=DataArrayDouble::New(); field->setArray(arr);
372 { arr->alloc(0,1); return; }
375 arr->alloc(nbnodes-1,1);
376 std::transform(_coords->begin()+1,_coords->end(),_coords->begin(),arr->getPointer(),std::minus<double>());
382 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> tmp=DataArrayDouble::New(); tmp->alloc(nbnodes-1,spaceDim);
383 std::transform(_coords->begin()+spaceDim,_coords->end(),_coords->begin(),tmp->getPointer(),std::minus<double>());
384 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> tmp2=tmp->magnitude(); field->setArray(tmp2);
389 * \param [in,out] f field feeded with good values.
390 * \sa MEDCouplingCurveLinearMesh::getMeasureField
392 void MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim2(bool isAbs, MEDCouplingFieldDouble *field) const throw(INTERP_KERNEL::Exception)
394 int nbcells=getNumberOfCells();
395 int spaceDim=getSpaceDimension();
396 if(spaceDim!=2 && spaceDim!=3)
397 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim2 : with meshDim 2 only space dimension 2 and 3 are possible !");
398 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=DataArrayDouble::New(); field->setArray(arr);
399 arr->alloc(nbcells,1);
400 double *pt=arr->getPointer();
401 const double *coords=_coords->begin();
402 int nX=_structure[0]-1;
404 for(int i=0;i<nbcells;i++,pt++)
406 int cy=i/nX,cx=i-cy*nX;
407 conn[0]=cy*(nX+1)+cx; conn[1]=(cy+1)*(nX+1)+cx; conn[2]=(cy+1)*(nX+1)+1+cx; conn[3]=cy*(nX+1)+cx+1;
408 *pt=INTERP_KERNEL::computeVolSurfOfCell2<int,INTERP_KERNEL::ALL_C_MODE>(INTERP_KERNEL::NORM_QUAD4,conn,4,coords,spaceDim);
415 * \param [in,out] f field feeded with good values.
416 * \sa MEDCouplingCurveLinearMesh::getMeasureField
418 void MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim3(bool isAbs, MEDCouplingFieldDouble *field) const throw(INTERP_KERNEL::Exception)
420 int nbcells=getNumberOfCells();
421 int spaceDim=getSpaceDimension();
423 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim3 : with meshDim 3 only space dimension 3 is possible !");
424 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=DataArrayDouble::New(); field->setArray(arr);
425 arr->alloc(nbcells,1);
426 double *pt=arr->getPointer();
427 const double *coords=_coords->begin();
428 int nX=_structure[0]-1,nY=(_structure[0]-1)*(_structure[1]-1);
429 int nY1=_structure[0]*_structure[1];
431 for(int i=0;i<nbcells;i++,pt++)
435 int cx=(i-cz*nY)-nX*cy;
436 conn[0]=cz*nY1+cy*(nX+1)+cx; conn[1]=cz*nY1+(cy+1)*(nX+1)+cx; conn[2]=cz*nY1+(cy+1)*(nX+1)+1+cx; conn[3]=cz*nY1+cy*(nX+1)+cx+1;
437 conn[4]=(cz+1)*nY1+cy*(nX+1)+cx; conn[5]=(cz+1)*nY1+(cy+1)*(nX+1)+cx; conn[6]=(cz+1)*nY1+(cy+1)*(nX+1)+1+cx; conn[7]=(cz+1)*nY1+cy*(nX+1)+cx+1;
438 *pt=INTERP_KERNEL::computeVolSurfOfCell2<int,INTERP_KERNEL::ALL_C_MODE>(INTERP_KERNEL::NORM_HEXA8,conn,8,coords,3);
445 * not implemented yet !
447 MEDCouplingFieldDouble *MEDCouplingCurveLinearMesh::getMeasureFieldOnNode(bool isAbs) const
449 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getMeasureFieldOnNode : not implemented yet !");
452 MEDCouplingFieldDouble *MEDCouplingCurveLinearMesh::buildOrthogonalField() const
454 if(getMeshDimension()!=2)
455 throw INTERP_KERNEL::Exception("Expected a cmesh with meshDim == 2 !");
456 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
457 DataArrayDouble *array=DataArrayDouble::New();
458 int nbOfCells=getNumberOfCells();
459 array->alloc(nbOfCells,3);
460 double *vals=array->getPointer();
461 for(int i=0;i<nbOfCells;i++)
462 { vals[3*i]=0.; vals[3*i+1]=0.; vals[3*i+2]=1.; }
463 ret->setArray(array);
473 template<const int SPACEDIMM>
477 static const int MY_SPACEDIM=SPACEDIMM;
478 static const int MY_MESHDIM=8;
479 typedef int MyConnType;
480 static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
482 // useless, but for windows compilation ...
483 const double* getCoordinatesPtr() const { return 0; }
484 const int* getConnectivityPtr() const { return 0; }
485 const int* getConnectivityIndexPtr() const { return 0; }
486 INTERP_KERNEL::NormalizedCellType getTypeOfElement(int) const { return (INTERP_KERNEL::NormalizedCellType)0; }
493 int MEDCouplingCurveLinearMesh::getCellContainingPoint(const double *pos, double eps) const
496 int spaceDim=getSpaceDimension();
497 const double *coords=_coords->getConstPointer();
499 _coords->distanceToTuple(pos,pos+spaceDim,nodeId);
501 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getCellContainingPoint : internal problem 1 !");
503 int nbOfNodes=getNumberOfNodes();
505 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getCellContainingPoint : No cells in this !");
506 switch(getMeshDimension())
513 conn[0]=nodeId-1; conn[1]=nodeId;
514 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<1> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_SEG2,coords,conn,2,eps))
517 if(nodeId<nbOfNodes-1)
519 conn[0]=nodeId; conn[1]=nodeId+1;
520 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<1> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_SEG2,coords,conn,2,eps))
527 int ny=nodeId/_structure[0],nx=nodeId-ny*_structure[0];
530 conn[0]=nx-1+_structure[0]*(ny-1); conn[1]=nx-1+_structure[0]*ny; conn[2]=nx+_structure[0]*ny; conn[3]=nx+_structure[0]*(ny-1);
531 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<2> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_QUAD4,coords,conn,4,eps))
532 return nx-1+(ny-1)*_structure[0];
534 if(nx<_structure[0]-1 && ny>0)
536 conn[0]=nx+_structure[0]*(ny-1); conn[1]=nx+_structure[0]*ny; conn[2]=nx+1+_structure[0]*ny; conn[3]=nx+1+_structure[0]*(ny-1);
537 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<2> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_QUAD4,coords,conn,4,eps))
538 return nx+(ny-1)*_structure[0];
540 if(nx>0 && ny<_structure[1]-1)
542 conn[0]=nx-1+_structure[0]*ny; conn[1]=nx-1+_structure[0]*(ny+1); conn[2]=nx+_structure[0]*(ny+1); conn[3]=nx+_structure[0]*ny;
543 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<2> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_QUAD4,coords,conn,4,eps))
544 return nx-1+ny*_structure[0];
546 if(nx<_structure[0]-1 && ny<_structure[1]-1)
548 conn[0]=nx+_structure[0]*ny; conn[1]=nx+_structure[0]*(ny+1); conn[2]=nx+1+_structure[0]*(ny+1); conn[3]=nx+1+_structure[0]*ny;
549 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<2> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_QUAD4,coords,conn,4,eps))
550 return nx+ny*_structure[0];
557 int nY=_structure[0]*_structure[1];
558 int nz=nodeId/_structure[1]; int ny=(nodeId-nz*nY)/_structure[0]; int nx=(nodeId-nz*nY)-_structure[0]*ny;
559 if(nx>0 && ny>0 && nz>0)
561 conn[0]=nx-1+_structure[0]*(ny-1)+nY*(nz-1); conn[1]=nx-1+_structure[2]*ny+nY*(nz-1); conn[2]=nx+_structure[2]*ny+nY*(nz-1); conn[3]=nx+_structure[0]*(ny-1)+nY*(nz-1);
562 conn[4]=nx-1+_structure[0]*(ny-1)+nY*nz; conn[5]=nx-1+_structure[0]*ny+nY*nz; conn[6]=nx+_structure[0]*ny+nY*nz; conn[7]=nx+_structure[0]*(ny-1)+nY*nz;
563 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
564 return nx-1+(ny-1)*_structure[0]+(nz-1)*nY;
566 if(nx<_structure[0]-1 && ny>0 && nz>0)
568 conn[0]=nx+_structure[0]*(ny-1)+nY*(nz-1); conn[1]=nx+_structure[0]*ny+nY*(nz-1); conn[2]=nx+1+_structure[0]*ny+nY*(nz-1); conn[3]=nx+1+_structure[0]*(ny-1)+nY*(nz-1);
569 conn[4]=nx+_structure[0]*(ny-1)+nY*nz; conn[5]=nx+_structure[0]*ny+nY*nz; conn[6]=nx+1+_structure[0]*ny+nY*nz; conn[7]=nx+1+_structure[0]*(ny-1)+nY*nz;
570 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
571 return nx+(ny-1)*_structure[0]+(nz-1)*nY;
573 if(nx>0 && ny<_structure[1]-1 && nz>0)
575 conn[0]=nx-1+_structure[0]*ny+nY*(nz-1); conn[1]=nx-1+_structure[0]*(ny+1)+nY*(nz-1); conn[2]=nx+_structure[0]*(ny+1)+nY*(nz-1); conn[3]=nx+_structure[0]*ny+nY*(nz-1);
576 conn[4]=nx-1+_structure[0]*ny+nY*nz; conn[5]=nx-1+_structure[0]*(ny+1)+nY*nz; conn[6]=nx+_structure[0]*(ny+1)+nY*nz; conn[7]=nx+_structure[0]*ny+nY*nz;
577 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
578 return nx-1+ny*_structure[0]+(nz-1)*nY;
580 if(nx<_structure[0]-1 && ny<_structure[1]-1 && nz>0)
582 conn[0]=nx+_structure[0]*ny+nY*(nz-1); conn[1]=nx+_structure[0]*(ny+1)+nY*(nz-1); conn[2]=nx+1+_structure[0]*(ny+1)+nY*(nz-1); conn[3]=nx+1+_structure[0]*ny+nY*(nz-1);
583 conn[4]=nx+_structure[0]*ny+nY*nz; conn[5]=nx+_structure[0]*(ny+1)+nY*nz; conn[6]=nx+1+_structure[0]*(ny+1)+nY*nz; conn[7]=nx+1+_structure[0]*ny+nY*nz;
584 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
585 return nx+ny*_structure[0]+(nz-1)*nY;
587 if(nx>0 && ny>0 && nz<_structure[2]-1)
589 conn[0]=nx-1+_structure[0]*(ny-1)+nY*(nz-1); conn[1]=nx-1+_structure[2]*ny+nY*(nz-1); conn[2]=nx+_structure[2]*ny+nY*(nz-1); conn[3]=nx+_structure[0]*(ny-1)+nY*(nz-1);
590 conn[4]=nx-1+_structure[0]*(ny-1)+nY*nz; conn[5]=nx-1+_structure[0]*ny+nY*nz; conn[6]=nx+_structure[0]*ny+nY*nz; conn[7]=nx+_structure[0]*(ny-1)+nY*nz;
591 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
592 return nx-1+(ny-1)*_structure[0]+nz*nY;
594 if(nx<_structure[0]-1 && ny>0 && nz<_structure[2]-1)
596 conn[0]=nx+_structure[0]*(ny-1)+nY*(nz-1); conn[1]=nx+_structure[0]*ny+nY*(nz-1); conn[2]=nx+1+_structure[0]*ny+nY*(nz-1); conn[3]=nx+1+_structure[0]*(ny-1)+nY*(nz-1);
597 conn[4]=nx+_structure[0]*(ny-1)+nY*nz; conn[5]=nx+_structure[0]*ny+nY*nz; conn[6]=nx+1+_structure[0]*ny+nY*nz; conn[7]=nx+1+_structure[0]*(ny-1)+nY*nz;
598 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
599 return nx+(ny-1)*_structure[0]+nz*nY;
601 if(nx>0 && ny<_structure[1]-1 && nz<_structure[2]-1)
603 conn[0]=nx-1+_structure[0]*ny+nY*(nz-1); conn[1]=nx-1+_structure[0]*(ny+1)+nY*(nz-1); conn[2]=nx+_structure[0]*(ny+1)+nY*(nz-1); conn[3]=nx+_structure[0]*ny+nY*(nz-1);
604 conn[4]=nx-1+_structure[0]*ny+nY*nz; conn[5]=nx-1+_structure[0]*(ny+1)+nY*nz; conn[6]=nx+_structure[0]*(ny+1)+nY*nz; conn[7]=nx+_structure[0]*ny+nY*nz;
605 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
606 return nx-1+ny*_structure[0]+nz*nY;
608 if(nx<_structure[0]-1 && ny<_structure[1]-1 && nz<_structure[2]-1)
610 conn[0]=nx+_structure[0]*ny+nY*(nz-1); conn[1]=nx+_structure[0]*(ny+1)+nY*(nz-1); conn[2]=nx+1+_structure[0]*(ny+1)+nY*(nz-1); conn[3]=nx+1+_structure[0]*ny+nY*(nz-1);
611 conn[4]=nx+_structure[0]*ny+nY*nz; conn[5]=nx+_structure[0]*(ny+1)+nY*nz; conn[6]=nx+1+_structure[0]*(ny+1)+nY*nz; conn[7]=nx+1+_structure[0]*ny+nY*nz;
612 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
613 return nx+ny*_structure[0]+nz*nY;
618 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getCellContainingPoint : mesh dimension managed are 1, 2 or 3 !");
622 void MEDCouplingCurveLinearMesh::rotate(const double *center, const double *vector, double angle)
624 if(!((DataArrayDouble *)_coords))
625 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::rotate : no coordinates set !");
626 int spaceDim=getSpaceDimension();
627 int nbNodes=_coords->getNumberOfTuples();
628 double *coords=_coords->getPointer();
630 MEDCouplingPointSet::Rotate3DAlg(center,vector,angle,nbNodes,coords);
632 MEDCouplingPointSet::Rotate2DAlg(center,angle,nbNodes,coords);
634 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::rotate : invalid space dim for rotation must be 2 or 3");
635 _coords->declareAsNew();
639 void MEDCouplingCurveLinearMesh::translate(const double *vector)
641 if(!((DataArrayDouble *)_coords))
642 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::translate : no coordinates set !");
643 double *coords=_coords->getPointer();
644 int nbNodes=getNumberOfNodes();
645 int dim=getSpaceDimension();
646 for(int i=0; i<nbNodes; i++)
647 for(int idim=0; idim<dim;idim++)
648 coords[i*dim+idim]+=vector[idim];
649 _coords->declareAsNew();
653 void MEDCouplingCurveLinearMesh::scale(const double *point, double factor)
655 if(!((DataArrayDouble *)_coords))
656 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::scale : no coordinates set !");
657 double *coords=_coords->getPointer();
658 int nbNodes=_coords->getNumberOfTuples();
659 int dim=_coords->getNumberOfComponents();
660 for(int i=0;i<nbNodes;i++)
662 std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
663 std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind2nd(std::multiplies<double>(),factor));
664 std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
666 _coords->declareAsNew();
670 MEDCouplingMesh *MEDCouplingCurveLinearMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
672 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::mergeMyselfWith : not available for CurveLinear Mesh !");
675 DataArrayDouble *MEDCouplingCurveLinearMesh::getCoordinatesAndOwner() const
677 DataArrayDouble *ret=const_cast<DataArrayDouble *>((const DataArrayDouble *)_coords);
683 DataArrayDouble *MEDCouplingCurveLinearMesh::getBarycenterAndOwner() const
686 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
687 int spaceDim=getSpaceDimension();
688 int meshDim=getMeshDimension();
689 int nbOfCells=getNumberOfCells();
690 ret->alloc(nbOfCells,spaceDim);
691 ret->copyStringInfoFrom(*getCoords());
695 { getBarycenterAndOwnerMeshDim3(ret); return ret.retn(); }
697 { getBarycenterAndOwnerMeshDim2(ret); return ret.retn(); }
699 { getBarycenterAndOwnerMeshDim1(ret); return ret.retn(); }
701 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getBarycenterAndOwner : mesh dimension must be in [1,2,3] !");
706 * \param [in,out] bary Barycenter array feeded with good values.
707 * \sa MEDCouplingCurveLinearMesh::getBarycenterAndOwner
709 void MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim3(DataArrayDouble *bary) const
711 int nbOfCells=getNumberOfCells();
712 double *ptToFill=bary->getPointer();
713 const double *coor=_coords->getConstPointer();
714 if(getSpaceDimension()!=3)
715 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim3 : with meshDim 3 only space dimension 3 is possible !");
716 int nX=_structure[0]-1,nY=(_structure[0]-1)*(_structure[1]-1);
717 int nY1=_structure[0]*_structure[1];
719 for(int i=0;i<nbOfCells;i++)
723 int cx=(i-cz*nY)-nX*cy;
724 conn[0]=cz*nY1+cy*(nX+1)+cx+1; conn[1]=cz*nY1+cy*(nX+1)+cx; conn[2]=cz*nY1+(cy+1)*(nX+1)+cx; conn[3]=cz*nY1+(cy+1)*(nX+1)+1+cx;
725 conn[4]=(cz+1)*nY1+cy*(nX+1)+cx+1; conn[5]=(cz+1)*nY1+cy*(nX+1)+cx; conn[6]=(cz+1)*nY1+(cy+1)*(nX+1)+cx; conn[7]=(cz+1)*nY1+(cy+1)*(nX+1)+1+cx;
726 INTERP_KERNEL::computeBarycenter2<int,INTERP_KERNEL::ALL_C_MODE>(INTERP_KERNEL::NORM_HEXA8,conn,8,coor,3,ptToFill);
732 * \param [in,out] bary Barycenter array feeded with good values.
733 * \sa MEDCouplingCurveLinearMesh::getBarycenterAndOwner
735 void MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim2(DataArrayDouble *bary) const
737 int nbcells=getNumberOfCells();
738 int spaceDim=getSpaceDimension();
739 double *ptToFill=bary->getPointer();
740 const double *coor=_coords->getConstPointer();
741 if(spaceDim!=2 && spaceDim!=3)
742 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim2 : with meshDim 2 only space dimension 2 and 3 are possible !");
743 int nX=_structure[0]-1;
745 for(int i=0;i<nbcells;i++)
747 int cy=i/nX,cx=i-cy*nX;
748 conn[0]=cy*(nX+1)+cx; conn[1]=(cy+1)*(nX+1)+cx; conn[2]=(cy+1)*(nX+1)+1+cx; conn[3]=cy*(nX+1)+cx+1;
749 INTERP_KERNEL::computeBarycenter2<int,INTERP_KERNEL::ALL_C_MODE>(INTERP_KERNEL::NORM_QUAD4,conn,4,coor,spaceDim,ptToFill);
755 * \param [in,out] bary Barycenter array feeded with good values.
756 * \sa MEDCouplingCurveLinearMesh::getBarycenterAndOwner
758 void MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim1(DataArrayDouble *bary) const
760 int spaceDim=getSpaceDimension();
761 std::transform(_coords->begin()+spaceDim,_coords->end(),_coords->begin(),bary->getPointer(),std::plus<double>());
762 std::transform(bary->begin(),bary->end(),bary->getPointer(),std::bind2nd(std::multiplies<double>(),0.5));
765 void MEDCouplingCurveLinearMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
767 throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CurveLinear Mesh !");
770 void MEDCouplingCurveLinearMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
773 double time=getTime(it,order);
776 littleStrings.clear();
777 littleStrings.push_back(getName());
778 littleStrings.push_back(getDescription());
779 littleStrings.push_back(getTimeUnit());
781 std::vector<std::string> littleStrings2;
782 if((const DataArrayDouble *)_coords)
783 _coords->getTinySerializationStrInformation(littleStrings2);
784 littleStrings.insert(littleStrings.end(),littleStrings2.begin(),littleStrings2.end());
786 tinyInfo.push_back(it);
787 tinyInfo.push_back(order);
788 tinyInfo.push_back((int)_structure.size());
789 for(std::vector<int>::const_iterator itt=_structure.begin();itt!=_structure.end();itt++)
790 tinyInfo.push_back(*itt);
791 std::vector<int> tinyInfo2;
792 if((const DataArrayDouble *)_coords)
793 _coords->getTinySerializationIntInformation(tinyInfo2);
794 tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
796 tinyInfoD.push_back(time);
799 void MEDCouplingCurveLinearMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
801 a1->alloc(tinyInfo[2],1);
802 std::vector<int> tinyInfo2(tinyInfo.begin()+3+tinyInfo[2],tinyInfo.end());
803 a2->resizeForUnserialization(tinyInfo2);
806 void MEDCouplingCurveLinearMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
808 a1=DataArrayInt::New();
809 a1->alloc((int)_structure.size(),1);
810 int *ptr=a1->getPointer();
811 for(std::vector<int>::const_iterator it=_structure.begin();it!=_structure.end();it++,ptr++)
814 if((const DataArrayDouble *)_coords)
815 if(_coords->isAllocated())
816 sz=_coords->getNbOfElems();
817 a2=DataArrayDouble::New();
819 if(sz!=0 && (const DataArrayDouble *)_coords)
820 std::copy(_coords->begin(),_coords->end(),a2->getPointer());
823 void MEDCouplingCurveLinearMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
824 const std::vector<std::string>& littleStrings)
826 setName(littleStrings[0].c_str());
827 setDescription(littleStrings[1].c_str());
828 setTimeUnit(littleStrings[2].c_str());
829 setTime(tinyInfoD[0],tinyInfo[0],tinyInfo[1]);
831 _structure.resize(sz);
832 for(int i=0;i<sz;i++)
833 _structure[i]=tinyInfo[3+i];
834 if((int)tinyInfo.size()>sz+3)
836 _coords=DataArrayDouble::New();
837 std::vector<int> tinyInfo2(tinyInfo.begin()+3+sz,tinyInfo.end());
838 _coords->resizeForUnserialization(tinyInfo2);
839 std::copy(a2->begin(),a2->end(),_coords->getPointer());
840 std::vector<std::string> littleStrings2(littleStrings.begin()+3,littleStrings.end());
841 _coords->finishUnserialization(tinyInfo2,littleStrings2);
845 void MEDCouplingCurveLinearMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception)
847 std::ostringstream extent;
848 int meshDim=(int)_structure.size();
849 if(meshDim<=0 || meshDim>3)
850 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::writeVTKLL : meshDim invalid ! must be in [1,2,3] !");
852 { int val=i<meshDim?_structure[i]-1:0; extent << "0 " << val << " "; }
853 ofs << " <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\">\n";
854 ofs << " <Piece Extent=\"" << extent.str() << "\">\n";
855 ofs << " <PointData>\n" << pointData << std::endl;
856 ofs << " </PointData>\n";
857 ofs << " <CellData>\n" << cellData << std::endl;
858 ofs << " </CellData>\n";
859 ofs << " <Points>\n";
860 if(getSpaceDimension()==3)
861 _coords->writeVTK(ofs,8,"Points");
864 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=_coords->changeNbOfComponents(3,0.);
865 coo->writeVTK(ofs,8,"Points");
867 ofs << " </Points>\n";
868 ofs << " </Piece>\n";
869 ofs << " </" << getVTKDataSetType() << ">\n";
872 std::string MEDCouplingCurveLinearMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception)
874 return std::string("StructuredGrid");