1 // Copyright (C) 2007-2023 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (EDF R&D)
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 MEDCoupling;
36 MEDCouplingCurveLinearMesh::MEDCouplingCurveLinearMesh():_coords(0),_structure(0)
40 MEDCouplingCurveLinearMesh::MEDCouplingCurveLinearMesh(const MEDCouplingCurveLinearMesh& other, bool deepCpy):MEDCouplingStructuredMesh(other,deepCpy),_structure(other._structure)
44 if((const DataArrayDouble *)other._coords)
45 _coords=other._coords->deepCopy();
48 _coords=other._coords;
51 MEDCouplingCurveLinearMesh::~MEDCouplingCurveLinearMesh()
55 MEDCouplingCurveLinearMesh *MEDCouplingCurveLinearMesh::New()
57 return new MEDCouplingCurveLinearMesh;
60 MEDCouplingCurveLinearMesh *MEDCouplingCurveLinearMesh::New(const std::string& meshName)
62 MEDCouplingCurveLinearMesh *ret=new MEDCouplingCurveLinearMesh;
63 ret->setName(meshName);
67 MEDCouplingCurveLinearMesh *MEDCouplingCurveLinearMesh::deepCopy() 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::getHeapMemorySizeWithoutChildren() const
85 std::size_t ret(MEDCouplingStructuredMesh::getHeapMemorySizeWithoutChildren());
86 ret+=_structure.capacity()*sizeof(mcIdType);
90 std::vector<const BigMemoryObject *> MEDCouplingCurveLinearMesh::getDirectChildrenWithNull() const
92 std::vector<const BigMemoryObject *> ret;
93 ret.push_back((const DataArrayDouble *)_coords);
98 * This method copyies all tiny strings from other (name and components name).
99 * @throw if other and this have not same mesh type.
101 void MEDCouplingCurveLinearMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
103 const MEDCouplingCurveLinearMesh *otherC=dynamic_cast<const MEDCouplingCurveLinearMesh *>(other);
105 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::copyTinyStringsFrom : meshes have not same type !");
106 MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
107 if((DataArrayDouble *)_coords && (const DataArrayDouble *)otherC->_coords)
108 _coords->copyStringInfoFrom(*otherC->_coords);
111 bool MEDCouplingCurveLinearMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
114 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::isEqualIfNotWhy : input other pointer is null !");
115 const MEDCouplingCurveLinearMesh *otherC=dynamic_cast<const MEDCouplingCurveLinearMesh *>(other);
118 reason="mesh given in input is not castable in MEDCouplingCurveLinearMesh !";
121 if(!MEDCouplingStructuredMesh::isEqualIfNotWhy(other,prec,reason))
123 std::ostringstream oss; oss.precision(15);
124 if(((const DataArrayDouble *)_coords && ((const DataArrayDouble *)otherC->_coords)==0) || (((const DataArrayDouble *)_coords)==0 && (const DataArrayDouble *)otherC->_coords))
126 oss << "Only one CurveLinearMesh between the two this and other has its coordinates defined !";
130 if((const DataArrayDouble *)_coords)
132 if(!_coords->isEqualIfNotWhy(*(otherC->_coords),prec,reason))
134 oss << "Coordinates DataArrayDouble of differ :";
135 reason.insert(0,oss.str());
138 if(_structure!=otherC->_structure)
139 { reason="CurveLinearMesh structures differ !"; return false; }
144 bool MEDCouplingCurveLinearMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
146 const MEDCouplingCurveLinearMesh *otherC=dynamic_cast<const MEDCouplingCurveLinearMesh *>(other);
149 if(((const DataArrayDouble *)_coords && ((const DataArrayDouble *)otherC->_coords)==0) || (((const DataArrayDouble *)_coords)==0 && (const DataArrayDouble *)otherC->_coords))
151 if((const DataArrayDouble *)_coords)
153 if(!_coords->isEqualWithoutConsideringStr(*(otherC->_coords),prec))
155 if(_structure!=otherC->_structure)
161 void MEDCouplingCurveLinearMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
162 DataArrayIdType *&cellCor, DataArrayIdType *&nodeCor) const
164 if(!isEqualWithoutConsideringStr(other,prec))
165 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkDeepEquivalWith : Meshes are not the same !");
169 * Nothing is done here (except to check that the other is a MEDCoupling::MEDCouplingCurveLinearMesh instance too).
170 * The user intend that the nodes are the same, so by construction of MEDCoupling::MEDCouplingCurveLinearMesh, \a this and \a other are the same !
172 void MEDCouplingCurveLinearMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
173 DataArrayIdType *&cellCor) const
175 if(!isEqualWithoutConsideringStr(other,prec))
176 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkDeepEquivalOnSameNodesWith : Meshes are not the same !");
179 void MEDCouplingCurveLinearMesh::checkConsistencyLight() const
181 std::size_t sz=_structure.size(),i=0;
182 mcIdType nbOfNodes=1;
184 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkConsistencyLight : structure should have a lgth of size 1 at least !");
185 for(std::vector<mcIdType>::const_iterator it=_structure.begin();it!=_structure.end();it++,i++)
188 { std::ostringstream oss; oss << "MEDCouplingCurveLinearMesh::checkConsistencyLight : At pos #" << i << " of structure value is " << *it << "should be >= 1 !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
191 if(!((const DataArrayDouble *)_coords))
192 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkConsistencyLight : the array is not set !");
193 if(!_coords->isAllocated())
194 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkConsistencyLight : the array is not allocated !");
195 if(_coords->getNumberOfComponents()<1)
196 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::checkConsistencyLight : the array should have >= 1 components !");
197 if(_coords->getNumberOfTuples()!=nbOfNodes)
199 std::ostringstream oss; oss << "MEDCouplingCurveLinearMesh::checkConsistencyLight : structure said that number of nodes should be equal to " << nbOfNodes << " but number of tuples in array is equal to " << _coords->getNumberOfTuples() << " !";
200 throw INTERP_KERNEL::Exception(oss.str().c_str());
204 void MEDCouplingCurveLinearMesh::checkConsistency(double eps) const
206 checkConsistencyLight();
209 mcIdType MEDCouplingCurveLinearMesh::getNumberOfCells() const
211 checkConsistencyLight();
212 return MEDCouplingStructuredMesh::getNumberOfCells();
215 mcIdType MEDCouplingCurveLinearMesh::getNumberOfNodes() const
217 checkConsistencyLight();
218 return MEDCouplingStructuredMesh::getNumberOfNodes();
221 void MEDCouplingCurveLinearMesh::getNodeGridStructure(mcIdType *res) const
223 std::copy(_structure.begin(),_structure.end(),res);
227 * MEDCouplingCurveLinearMesh has the property to define 2 space dimensions. One coming from its coordinates. The other coming from the node structure.
228 * Normally they should be equal ! This method returns the space dimension from coordinates. If the other one is requested call getSpaceDimensionOnNodeStruct.
230 * \sa MEDCouplingStructuredMesh::getSpaceDimensionOnNodeStruct
232 int MEDCouplingCurveLinearMesh::getSpaceDimension() const
234 if(!((const DataArrayDouble *)_coords))
235 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getSpaceDimension : no array set ! impossible to deduce a space dimension !");
236 return int(_coords->getNumberOfComponents());
239 void MEDCouplingCurveLinearMesh::getCoordinatesOfNode(mcIdType nodeId, std::vector<double>& coo) const
241 if(!((const DataArrayDouble *)_coords))
242 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getCoordinatesOfNode : Coordinates not set !");
243 std::size_t nbOfCompo=_coords->getNumberOfComponents();
244 if(nodeId>=0 && nodeId<_coords->getNumberOfTuples())
245 coo.insert(coo.end(),_coords->begin()+nodeId*nbOfCompo,_coords->begin()+(nodeId+1)*nbOfCompo);
247 { std::ostringstream oss; oss << "MEDCouplingCurveLinearMesh::getCoordinatesOfNode : nodeId has to be in [0," << _coords->getNumberOfTuples() << ") !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
250 std::string MEDCouplingCurveLinearMesh::simpleRepr() const
252 std::ostringstream ret;
253 ret << "Curve linear mesh with name : \"" << getName() << "\"\n";
254 ret << "Description of mesh : \"" << getDescription() << "\"\n";
256 double tt=getTime(tmpp1,tmpp2);
257 ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
258 ret << "Iteration : " << tmpp1 << " Order : " << tmpp2 << "\n";
259 ret << "The nodal structure of curve linear mesh is : [";
260 std::copy(_structure.begin(),_structure.end(),std::ostream_iterator<int>(ret,",")); ret << "]\n";
261 ret << "The coords array is this : ";
262 if((const DataArrayDouble *)_coords)
263 _coords->reprZipWithoutNameStream(ret);
265 ret << "no array specified !";
269 std::string MEDCouplingCurveLinearMesh::advancedRepr() const
274 const DataArrayDouble *MEDCouplingCurveLinearMesh::getDirectAccessOfCoordsArrIfInStructure() const
279 DataArrayDouble *MEDCouplingCurveLinearMesh::getCoords()
284 const DataArrayDouble *MEDCouplingCurveLinearMesh::getCoords() const
289 void MEDCouplingCurveLinearMesh::setCoords(const DataArrayDouble *coords)
291 if(coords!=(const DataArrayDouble *)_coords)
293 _coords=const_cast<DataArrayDouble *>(coords);
300 void MEDCouplingCurveLinearMesh::setNodeGridStructure(const mcIdType *gridStructBg, const mcIdType *gridStructEnd)
302 std::size_t sz=std::distance(gridStructBg,gridStructEnd);
305 _structure.resize(0);
306 _structure.insert(_structure.end(),gridStructBg,gridStructEnd);
310 std::ostringstream oss; oss << "MEDCouplingCurveLinearMesh::setNodeGridStructure : size of input nodal grid structure (" << sz << ") should be in 1, 2 or 3 !";
311 throw INTERP_KERNEL::Exception(oss.str().c_str());
315 std::vector<mcIdType> MEDCouplingCurveLinearMesh::getNodeGridStructure() const
320 MEDCouplingStructuredMesh *MEDCouplingCurveLinearMesh::buildStructuredSubPart(const std::vector< std::pair<mcIdType,mcIdType> >& cellPart) const
322 checkConsistencyLight();
323 int dim(getSpaceDimension());
324 std::vector<mcIdType> dims(getMeshDimension());
325 if(dim!=ToIdType(cellPart.size()))
327 std::ostringstream oss; oss << "MEDCouplingCurveLinearMesh::buildStructuredSubPart : the space dimension is " << dim << " and cell part size is " << cellPart.size() << " !";
328 throw INTERP_KERNEL::Exception(oss.str().c_str());
330 std::vector< std::pair<mcIdType,mcIdType> > nodePartFormat(cellPart);
331 for(std::vector< std::pair<mcIdType,mcIdType> >::iterator it=nodePartFormat.begin();it!=nodePartFormat.end();it++)
333 MCAuto<DataArrayIdType> tmp1(BuildExplicitIdsFrom(getNodeGridStructure(),nodePartFormat));
334 MCAuto<MEDCouplingCurveLinearMesh> ret(dynamic_cast<MEDCouplingCurveLinearMesh *>(deepCopy()));
335 const DataArrayDouble *coo(ret->getCoords());
338 MCAuto<DataArrayDouble> coo2(coo->selectByTupleIdSafe(tmp1->begin(),tmp1->end()));
339 ret->setCoords(coo2);
341 for(int i=0;i<dim;i++)
343 dims[i]=cellPart[i].second-cellPart[i].first+1;
345 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::buildStructuredSubPart : invalid input cellPart !");
347 ret->setNodeGridStructure(&dims[0],&dims[0]+dims.size());
351 void MEDCouplingCurveLinearMesh::getBoundingBox(double *bbox) const
353 if(!((const DataArrayDouble *)_coords))
354 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getBoundingBox : Coordinates not set !");
355 _coords->getMinMaxPerComponent(bbox);
358 MEDCouplingFieldDouble *MEDCouplingCurveLinearMesh::getMeasureField(bool isAbs) const
360 checkConsistencyLight();
361 int meshDim=getMeshDimension();
362 std::string name="MeasureOfMesh_"; name+=getName();
363 MCAuto<MEDCouplingFieldDouble> field=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
364 field->setName(name); field->setMesh(const_cast<MEDCouplingCurveLinearMesh *>(this)); field->synchronizeTimeWithMesh();
368 { getMeasureFieldMeshDim3(isAbs,field); return field.retn(); }
370 { getMeasureFieldMeshDim2(isAbs,field); return field.retn(); }
372 { getMeasureFieldMeshDim1(isAbs,field); return field.retn(); }
374 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getMeasureField : mesh dimension must be in [1,2,3] !");
379 * \param [in] isAbs whether to compute signed or absolute values
380 * \param [in,out] field field fed with good values.
381 * \sa MEDCouplingCurveLinearMesh::getMeasureField
383 void MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim1(bool isAbs, MEDCouplingFieldDouble *field) const
385 mcIdType nbnodes=getNumberOfNodes();
386 int spaceDim=getSpaceDimension();
387 MCAuto<DataArrayDouble> arr=DataArrayDouble::New(); field->setArray(arr);
389 { arr->alloc(0,1); return; }
392 arr->alloc(nbnodes-1,1);
393 std::transform(_coords->begin()+1,_coords->end(),_coords->begin(),arr->getPointer(),std::minus<double>());
399 MCAuto<DataArrayDouble> tmp=DataArrayDouble::New(); tmp->alloc(nbnodes-1,spaceDim);
400 std::transform(_coords->begin()+spaceDim,_coords->end(),_coords->begin(),tmp->getPointer(),std::minus<double>());
401 MCAuto<DataArrayDouble> tmp2=tmp->magnitude(); field->setArray(tmp2);
406 * \param [in] isAbs whether to compute signed or absolute values
407 * \param [in,out] field field fed with good values.
408 * \sa MEDCouplingCurveLinearMesh::getMeasureField
410 void MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim2(bool isAbs, MEDCouplingFieldDouble *field) const
412 mcIdType nbcells=getNumberOfCells();
413 int spaceDim=getSpaceDimension();
414 if(spaceDim!=2 && spaceDim!=3)
415 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim2 : with meshDim 2 only space dimension 2 and 3 are possible !");
416 MCAuto<DataArrayDouble> arr=DataArrayDouble::New(); field->setArray(arr);
417 arr->alloc(nbcells,1);
418 double *pt=arr->getPointer();
419 const double *coords=_coords->begin();
420 mcIdType nX=_structure[0]-1;
422 for(mcIdType i=0;i<nbcells;i++,pt++)
424 mcIdType cy=i/nX,cx=i-cy*nX;
425 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;
426 *pt=INTERP_KERNEL::computeVolSurfOfCell2<mcIdType,INTERP_KERNEL::ALL_C_MODE>(INTERP_KERNEL::NORM_QUAD4,conn,4,coords,spaceDim);
433 * \param [in] isAbs whether to compute signed or absolute values
434 * \param [in,out] field field fed with good values.
435 * \sa MEDCouplingCurveLinearMesh::getMeasureField
437 void MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim3(bool isAbs, MEDCouplingFieldDouble *field) const
439 mcIdType nbcells=getNumberOfCells();
440 int spaceDim=getSpaceDimension();
442 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getMeasureFieldMeshDim3 : with meshDim 3 only space dimension 3 is possible !");
443 MCAuto<DataArrayDouble> arr=DataArrayDouble::New(); field->setArray(arr);
444 arr->alloc(nbcells,1);
445 double *pt=arr->getPointer();
446 const double *coords=_coords->begin();
447 mcIdType nX=_structure[0]-1,nY=(_structure[0]-1)*(_structure[1]-1);
448 mcIdType nY1=_structure[0]*_structure[1];
450 for(mcIdType i=0;i<nbcells;i++,pt++)
453 mcIdType cy=(i-cz*nY)/nX;
454 mcIdType cx=(i-cz*nY)-nX*cy;
455 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;
456 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;
457 *pt=INTERP_KERNEL::computeVolSurfOfCell2<mcIdType,INTERP_KERNEL::ALL_C_MODE>(INTERP_KERNEL::NORM_HEXA8,conn,8,coords,3);
464 * not implemented yet !
466 MEDCouplingFieldDouble *MEDCouplingCurveLinearMesh::getMeasureFieldOnNode(bool isAbs) const
468 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getMeasureFieldOnNode : not implemented yet !");
471 MEDCouplingFieldDouble *MEDCouplingCurveLinearMesh::buildOrthogonalField() const
473 if(getMeshDimension()!=2)
474 throw INTERP_KERNEL::Exception("Expected a cmesh with meshDim == 2 !");
475 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
476 DataArrayDouble *array=DataArrayDouble::New();
477 mcIdType nbOfCells=getNumberOfCells();
478 array->alloc(nbOfCells,3);
479 double *vals=array->getPointer();
480 for(mcIdType i=0;i<nbOfCells;i++)
481 { vals[3*i]=0.; vals[3*i+1]=0.; vals[3*i+2]=1.; }
482 ret->setArray(array);
490 namespace MEDCoupling
492 template<const int SPACEDIMM>
496 static const int MY_SPACEDIM=SPACEDIMM;
497 static const int MY_MESHDIM=8;
498 typedef mcIdType MyConnType;
499 static const INTERP_KERNEL::NumberingPolicy My_numPol=INTERP_KERNEL::ALL_C_MODE;
501 // useless, but for windows compilation ...
502 const double* getCoordinatesPtr() const { return 0; }
503 const mcIdType* getConnectivityPtr() const { return 0; }
504 const mcIdType* getConnectivityIndexPtr() const { return 0; }
505 INTERP_KERNEL::NormalizedCellType getTypeOfElement(mcIdType) const { return (INTERP_KERNEL::NormalizedCellType)0; }
512 mcIdType MEDCouplingCurveLinearMesh::getCellContainingPoint(const double *pos, double eps) const
514 checkConsistencyLight();
515 int spaceDim=getSpaceDimension();
516 const double *coords=_coords->getConstPointer();
518 _coords->distanceToTuple(pos,pos+spaceDim,nodeId);
520 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getCellContainingPoint : internal problem 1 !");
522 mcIdType nbOfNodes=getNumberOfNodes();
524 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getCellContainingPoint : No cells in this !");
525 switch(getMeshDimension())
532 conn[0]=nodeId-1; conn[1]=nodeId;
533 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<1> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_SEG2,coords,conn,2,eps))
536 if(nodeId<nbOfNodes-1)
538 conn[0]=nodeId; conn[1]=nodeId+1;
539 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<1> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_SEG2,coords,conn,2,eps))
547 mcIdType ny=nodeId/_structure[0],nx=nodeId-ny*_structure[0];
550 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);
551 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<2> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_QUAD4,coords,conn,4,eps))
552 return nx-1+(ny-1)*_structure[0];
554 if(nx<_structure[0]-1 && ny>0)
556 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);
557 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<2> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_QUAD4,coords,conn,4,eps))
558 return nx+(ny-1)*_structure[0];
560 if(nx>0 && ny<_structure[1]-1)
562 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;
563 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<2> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_QUAD4,coords,conn,4,eps))
564 return nx-1+ny*_structure[0];
566 if(nx<_structure[0]-1 && ny<_structure[1]-1)
568 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;
569 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<2> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_QUAD4,coords,conn,4,eps))
570 return nx+ny*_structure[0];
578 mcIdType nY=_structure[0]*_structure[1];
579 mcIdType nz=nodeId/_structure[1]; mcIdType ny=(nodeId-nz*nY)/_structure[0]; mcIdType nx=(nodeId-nz*nY)-_structure[0]*ny;
580 if(nx>0 && ny>0 && nz>0)
582 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);
583 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;
584 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
585 return nx-1+(ny-1)*_structure[0]+(nz-1)*nY;
587 if(nx<_structure[0]-1 && ny>0 && nz>0)
589 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);
590 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;
591 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
592 return nx+(ny-1)*_structure[0]+(nz-1)*nY;
594 if(nx>0 && ny<_structure[1]-1 && nz>0)
596 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);
597 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;
598 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
599 return nx-1+ny*_structure[0]+(nz-1)*nY;
601 if(nx<_structure[0]-1 && ny<_structure[1]-1 && nz>0)
603 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);
604 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;
605 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
606 return nx+ny*_structure[0]+(nz-1)*nY;
608 if(nx>0 && ny>0 && nz<_structure[2]-1)
610 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);
611 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;
612 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
613 return nx-1+(ny-1)*_structure[0]+nz*nY;
615 if(nx<_structure[0]-1 && ny>0 && nz<_structure[2]-1)
617 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);
618 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;
619 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
620 return nx+(ny-1)*_structure[0]+nz*nY;
622 if(nx>0 && ny<_structure[1]-1 && nz<_structure[2]-1)
624 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);
625 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;
626 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
627 return nx-1+ny*_structure[0]+nz*nY;
629 if(nx<_structure[0]-1 && ny<_structure[1]-1 && nz<_structure[2]-1)
631 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);
632 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;
633 if(INTERP_KERNEL::PointLocatorAlgos<DummyClsMCL<3> >::isElementContainsPoint(pos,INTERP_KERNEL::NORM_HEXA8,coords,conn,8,eps))
634 return nx+ny*_structure[0]+nz*nY;
640 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getCellContainingPoint : mesh dimension managed are 1, 2 or 3 !");
645 void MEDCouplingCurveLinearMesh::getCellsContainingPoint(const double *pos, double eps, std::vector<mcIdType>& elts) const
647 mcIdType ret(getCellContainingPoint(pos,eps));
651 void MEDCouplingCurveLinearMesh::rotate(const double *center, const double *vector, double angle)
653 if(!((DataArrayDouble *)_coords))
654 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::rotate : no coordinates set !");
655 int spaceDim=getSpaceDimension();
656 mcIdType nbNodes(_coords->getNumberOfTuples());
657 double *coords=_coords->getPointer();
659 DataArrayDouble::Rotate3DAlg(center,vector,angle,nbNodes,coords,coords);
661 DataArrayDouble::Rotate2DAlg(center,angle,nbNodes,coords,coords);
663 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::rotate : invalid space dim for rotation must be 2 or 3");
664 _coords->declareAsNew();
668 void MEDCouplingCurveLinearMesh::translate(const double *vector)
671 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::translate : NULL input point !");
672 if(!((DataArrayDouble *)_coords))
673 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::translate : no coordinates set !");
674 double *coords=_coords->getPointer();
675 mcIdType nbNodes=getNumberOfNodes();
676 int dim=getSpaceDimension();
677 for(mcIdType i=0; i<nbNodes; i++)
678 for(int idim=0; idim<dim;idim++)
679 coords[i*dim+idim]+=vector[idim];
680 _coords->declareAsNew();
684 void MEDCouplingCurveLinearMesh::scale(const double *point, double factor)
687 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::scale : NULL input point !");
688 if(!((DataArrayDouble *)_coords))
689 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::scale : no coordinates set !");
690 double *coords=_coords->getPointer();
691 mcIdType nbNodes(_coords->getNumberOfTuples());
692 std::size_t dim(_coords->getNumberOfComponents());
693 for(mcIdType i=0;i<nbNodes;i++)
695 std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::minus<double>());
696 std::transform(coords+i*dim,coords+(i+1)*dim,coords+i*dim,std::bind(std::multiplies<double>(),std::placeholders::_1,factor));
697 std::transform(coords+i*dim,coords+(i+1)*dim,point,coords+i*dim,std::plus<double>());
699 _coords->declareAsNew();
703 MEDCouplingMesh *MEDCouplingCurveLinearMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
705 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::mergeMyselfWith : not available for CurveLinear Mesh !");
708 DataArrayDouble *MEDCouplingCurveLinearMesh::getCoordinatesAndOwner() const
710 DataArrayDouble *ret=const_cast<DataArrayDouble *>((const DataArrayDouble *)_coords);
716 DataArrayDouble *MEDCouplingCurveLinearMesh::computeCellCenterOfMass() const
718 checkConsistencyLight();
719 MCAuto<DataArrayDouble> ret=DataArrayDouble::New();
720 int spaceDim=getSpaceDimension();
721 int meshDim=getMeshDimension();
722 mcIdType nbOfCells=getNumberOfCells();
723 ret->alloc(nbOfCells,spaceDim);
724 ret->copyStringInfoFrom(*getCoords());
728 { getBarycenterAndOwnerMeshDim3(ret); return ret.retn(); }
730 { getBarycenterAndOwnerMeshDim2(ret); return ret.retn(); }
732 { getBarycenterAndOwnerMeshDim1(ret); return ret.retn(); }
734 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::computeCellCenterOfMass : mesh dimension must be in [1,2,3] !");
738 DataArrayDouble *MEDCouplingCurveLinearMesh::computeIsoBarycenterOfNodesPerCell() const
740 return MEDCouplingCurveLinearMesh::computeCellCenterOfMass();
744 * \param [in,out] bary Barycenter array fed with good values.
745 * \sa MEDCouplingCurveLinearMesh::computeCellCenterOfMass
747 void MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim3(DataArrayDouble *bary) const
749 mcIdType nbOfCells=getNumberOfCells();
750 double *ptToFill=bary->getPointer();
751 const double *coor=_coords->getConstPointer();
752 if(getSpaceDimension()!=3)
753 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim3 : with meshDim 3 only space dimension 3 is possible !");
754 mcIdType nX=_structure[0]-1,nY=(_structure[0]-1)*(_structure[1]-1);
755 mcIdType nY1=_structure[0]*_structure[1];
757 for(mcIdType i=0;i<nbOfCells;i++)
760 mcIdType cy=(i-cz*nY)/nX;
761 mcIdType cx=(i-cz*nY)-nX*cy;
762 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;
763 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;
764 INTERP_KERNEL::computeBarycenter2<mcIdType,INTERP_KERNEL::ALL_C_MODE>(INTERP_KERNEL::NORM_HEXA8,conn,8,coor,3,ptToFill);
770 * \param [in,out] bary Barycenter array fed with good values.
771 * \sa MEDCouplingCurveLinearMesh::computeCellCenterOfMass
773 void MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim2(DataArrayDouble *bary) const
775 mcIdType nbcells=getNumberOfCells();
776 int spaceDim=getSpaceDimension();
777 double *ptToFill=bary->getPointer();
778 const double *coor=_coords->getConstPointer();
779 if(spaceDim!=2 && spaceDim!=3)
780 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim2 : with meshDim 2 only space dimension 2 and 3 are possible !");
781 mcIdType nX=_structure[0]-1;
783 for(mcIdType i=0;i<nbcells;i++)
785 mcIdType cy=i/nX,cx=i-cy*nX;
786 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;
787 INTERP_KERNEL::computeBarycenter2<mcIdType,INTERP_KERNEL::ALL_C_MODE>(INTERP_KERNEL::NORM_QUAD4,conn,4,coor,spaceDim,ptToFill);
793 * \param [in,out] bary Barycenter array fed with good values.
794 * \sa MEDCouplingCurveLinearMesh::computeCellCenterOfMass
796 void MEDCouplingCurveLinearMesh::getBarycenterAndOwnerMeshDim1(DataArrayDouble *bary) const
798 int spaceDim=getSpaceDimension();
799 std::transform(_coords->begin()+spaceDim,_coords->end(),_coords->begin(),bary->getPointer(),std::plus<double>());
800 std::transform(bary->begin(),bary->end(),bary->getPointer(),std::bind(std::multiplies<double>(),std::placeholders::_1,0.5));
803 void MEDCouplingCurveLinearMesh::renumberCells(const mcIdType *old2NewBg, bool check)
805 throw INTERP_KERNEL::Exception("Functionality of renumbering cell not available for CurveLinear Mesh !");
808 void MEDCouplingCurveLinearMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<mcIdType>& tinyInfo, std::vector<std::string>& littleStrings) const
811 double time=getTime(it,order);
814 littleStrings.clear();
815 littleStrings.push_back(getName());
816 littleStrings.push_back(getDescription());
817 littleStrings.push_back(getTimeUnit());
819 std::vector<std::string> littleStrings2;
820 if((const DataArrayDouble *)_coords)
821 _coords->getTinySerializationStrInformation(littleStrings2);
822 littleStrings.insert(littleStrings.end(),littleStrings2.begin(),littleStrings2.end());
824 tinyInfo.push_back(it);
825 tinyInfo.push_back(order);
826 tinyInfo.push_back(ToIdType(_structure.size()));
827 for(std::vector<mcIdType>::const_iterator itt=_structure.begin();itt!=_structure.end();itt++)
828 tinyInfo.push_back(*itt);
829 std::vector<mcIdType> tinyInfo2;
830 if((const DataArrayDouble *)_coords)
831 _coords->getTinySerializationIntInformation(tinyInfo2);
832 tinyInfo.insert(tinyInfo.end(),tinyInfo2.begin(),tinyInfo2.end());
834 tinyInfoD.push_back(time);
837 void MEDCouplingCurveLinearMesh::resizeForUnserialization(const std::vector<mcIdType>& tinyInfo, DataArrayIdType *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
839 a1->alloc(tinyInfo[2],1);
840 std::vector<mcIdType> tinyInfo2(tinyInfo.begin()+3+tinyInfo[2],tinyInfo.end());
841 a2->resizeForUnserialization(tinyInfo2);
844 void MEDCouplingCurveLinearMesh::serialize(DataArrayIdType *&a1, DataArrayDouble *&a2) const
846 a1=DataArrayIdType::New();
847 a1->alloc(_structure.size(),1);
848 mcIdType *ptr=a1->getPointer();
849 for(std::vector<mcIdType>::const_iterator it=_structure.begin();it!=_structure.end();it++,ptr++)
852 if((const DataArrayDouble *)_coords)
853 if(_coords->isAllocated())
854 sz=_coords->getNbOfElems();
855 a2=DataArrayDouble::New();
857 if(sz!=0 && (const DataArrayDouble *)_coords)
858 std::copy(_coords->begin(),_coords->end(),a2->getPointer());
861 void MEDCouplingCurveLinearMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<mcIdType>& tinyInfo, const DataArrayIdType *a1, DataArrayDouble *a2,
862 const std::vector<std::string>& littleStrings)
864 setName(littleStrings[0]);
865 setDescription(littleStrings[1]);
866 setTimeUnit(littleStrings[2]);
867 setTime(tinyInfoD[0],FromIdType<int>(tinyInfo[0]),FromIdType<int>(tinyInfo[1]));
868 mcIdType sz=tinyInfo[2];
869 _structure.resize(sz);
870 for(mcIdType i=0;i<sz;i++)
871 _structure[i]=tinyInfo[3+i];
872 if(ToIdType(tinyInfo.size())>sz+3)
874 _coords=DataArrayDouble::New();
875 std::vector<mcIdType> tinyInfo2(tinyInfo.begin()+3+sz,tinyInfo.end());
876 _coords->resizeForUnserialization(tinyInfo2);
877 std::copy(a2->begin(),a2->end(),_coords->getPointer());
878 std::vector<std::string> littleStrings2(littleStrings.begin()+3,littleStrings.end());
879 _coords->finishUnserialization(tinyInfo2,littleStrings2);
883 void MEDCouplingCurveLinearMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
885 std::ostringstream extent;
886 std::size_t meshDim=_structure.size();
887 if(meshDim==0 || meshDim>3)
888 throw INTERP_KERNEL::Exception("MEDCouplingCurveLinearMesh::writeVTKLL : meshDim invalid ! must be in [1,2,3] !");
889 for(std::size_t i=0;i<3;i++)
890 { mcIdType val=i<meshDim?_structure[i]-1:0; extent << "0 " << val << " "; }
891 ofs << " <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\">\n";
892 ofs << " <Piece Extent=\"" << extent.str() << "\">\n";
893 ofs << " <PointData>\n" << pointData << std::endl;
894 ofs << " </PointData>\n";
895 ofs << " <CellData>\n" << cellData << std::endl;
896 ofs << " </CellData>\n";
897 ofs << " <Points>\n";
898 if(getSpaceDimension()==3)
899 _coords->writeVTK(ofs,8,"Points",byteData);
902 MCAuto<DataArrayDouble> coo=_coords->changeNbOfComponents(3,0.);
903 coo->writeVTK(ofs,8,"Points",byteData);
905 ofs << " </Points>\n";
906 ofs << " </Piece>\n";
907 ofs << " </" << getVTKDataSetType() << ">\n";
910 void MEDCouplingCurveLinearMesh::reprQuickOverview(std::ostream& stream) const
912 stream << "MEDCouplingCurveLinearMesh C++ instance at " << this << ". Name : \"" << getName() << "\".";
913 stream << " Nodal structure : [";
914 std::size_t s_size=_structure.size();
915 for(std::size_t i=0;i<s_size;i++)
917 char tmp=(char)((int)('X')+i);
918 stream << " " << tmp << "=" << _structure[i];
923 const DataArrayDouble *coo(_coords);
925 { stream << std::endl << "No coordinates set !"; return ; }
926 if(!coo->isAllocated())
927 { stream << std::endl << "Coordinates set but not allocated !"; return ; }
928 std::size_t nbOfCompo(coo->getNumberOfComponents());
929 std::size_t nbOfCompoExp(-1);
932 nbOfCompoExp=getSpaceDimension();
934 catch(INTERP_KERNEL::Exception&)
937 if(nbOfCompo!=nbOfCompoExp)
938 { stream << std::endl << "Coordinates set and allocated but mismatch number of components !"; return ; }
939 stream << std::endl << "Coordinates ( number of tuples = " << coo->getNumberOfTuples() << " ) : ";
940 coo->reprQuickOverviewData(stream,200);
943 std::string MEDCouplingCurveLinearMesh::getVTKFileExtension() const
945 return std::string("vts");
948 std::string MEDCouplingCurveLinearMesh::getVTKDataSetType() const
950 return std::string("StructuredGrid");