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 "MEDCouplingCMesh.hxx"
22 #include "MEDCouplingMemArray.hxx"
23 #include "MEDCouplingFieldDouble.hxx"
30 using namespace ParaMEDMEM;
32 MEDCouplingCMesh::MEDCouplingCMesh():_x_array(0),_y_array(0),_z_array(0)
36 MEDCouplingCMesh::MEDCouplingCMesh(const MEDCouplingCMesh& other, bool deepCopy):MEDCouplingStructuredMesh(other,deepCopy)
41 _x_array=other._x_array->deepCpy();
45 _y_array=other._y_array->deepCpy();
49 _z_array=other._z_array->deepCpy();
55 _x_array=other._x_array;
58 _y_array=other._y_array;
61 _z_array=other._z_array;
67 MEDCouplingCMesh::~MEDCouplingCMesh()
77 MEDCouplingCMesh *MEDCouplingCMesh::New()
79 return new MEDCouplingCMesh;
82 MEDCouplingCMesh *MEDCouplingCMesh::New(const char *meshName)
84 MEDCouplingCMesh *ret=new MEDCouplingCMesh;
85 ret->setName(meshName);
89 MEDCouplingMesh *MEDCouplingCMesh::deepCpy() const
94 MEDCouplingCMesh *MEDCouplingCMesh::clone(bool recDeepCpy) const
96 return new MEDCouplingCMesh(*this,recDeepCpy);
99 void MEDCouplingCMesh::updateTime() const
102 updateTimeWith(*_x_array);
104 updateTimeWith(*_y_array);
106 updateTimeWith(*_z_array);
109 std::size_t MEDCouplingCMesh::getHeapMemorySize() const
112 std::set<DataArrayDouble *> s;
113 s.insert(_x_array); s.insert(_y_array); s.insert(_z_array);
115 for(std::set<DataArrayDouble *>::const_iterator it=s.begin();it!=s.end();it++)
117 ret+=(*it)->getHeapMemorySize();
118 return MEDCouplingStructuredMesh::getHeapMemorySize()+ret;
122 * This method copyies all tiny strings from other (name and components name).
123 * @throw if other and this have not same mesh type.
125 void MEDCouplingCMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
127 const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
129 throw INTERP_KERNEL::Exception("MEDCouplingCMesh::copyTinyStringsFrom : meshes have not same type !");
130 MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
131 if(_x_array && otherC->_x_array)
132 _x_array->copyStringInfoFrom(*otherC->_x_array);
133 if(_y_array && otherC->_y_array)
134 _y_array->copyStringInfoFrom(*otherC->_y_array);
135 if(_z_array && otherC->_z_array)
136 _z_array->copyStringInfoFrom(*otherC->_z_array);
139 bool MEDCouplingCMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
142 throw INTERP_KERNEL::Exception("MEDCouplingCMesh::isEqualIfNotWhy : input other pointer is null !");
143 const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
146 reason="mesh given in input is not castable in MEDCouplingCMesh !";
149 if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
151 const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
152 const DataArrayDouble *otherArr[3]={otherC->_x_array,otherC->_y_array,otherC->_z_array};
153 std::ostringstream oss; oss.precision(15);
156 if((thisArr[i]!=0 && otherArr[i]==0) || (thisArr[i]==0 && otherArr[i]!=0))
158 oss << "Only one CMesh between the two this and other has its coordinates of rank" << i << " defined !";
163 if(!thisArr[i]->isEqualIfNotWhy(*otherArr[i],prec,reason))
165 oss << "Coordinates DataArrayDouble of rank #" << i << " differ :";
166 reason.insert(0,oss.str());
173 bool MEDCouplingCMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
175 const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
178 const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
179 const DataArrayDouble *otherArr[3]={otherC->_x_array,otherC->_y_array,otherC->_z_array};
182 if((thisArr[i]!=0 && otherArr[i]==0) || (thisArr[i]==0 && otherArr[i]!=0))
185 if(!thisArr[i]->isEqualWithoutConsideringStr(*otherArr[i],prec))
191 void MEDCouplingCMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
192 DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
194 if(!isEqualWithoutConsideringStr(other,prec))
195 throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkDeepEquivalWith : Meshes are not the same !");
199 * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingCMesh instance too).
200 * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingCMesh, 'this' and 'other' are the same !
202 void MEDCouplingCMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
203 DataArrayInt *&cellCor) const throw(INTERP_KERNEL::Exception)
205 const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
207 throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkDeepEquivalOnSameNodesWith : other is NOT a cartesian mesh ! Impossible to check equivalence !");
210 void MEDCouplingCMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
212 const char msg0[]="Invalid ";
213 const char msg1[]=" array ! Must contain more than 1 element.";
214 const char msg2[]=" array ! Must be with only one component.";
217 if(_x_array->getNbOfElems()<2)
219 std::ostringstream os; os << msg0 << 'X' << msg1;
220 throw INTERP_KERNEL::Exception(os.str().c_str());
222 if(_x_array->getNumberOfComponents()!=1)
224 std::ostringstream os; os << msg0 << 'X' << msg2;
225 throw INTERP_KERNEL::Exception(os.str().c_str());
230 if(_y_array->getNbOfElems()<2)
232 std::ostringstream os; os << msg0 << 'Y' << msg1;
233 throw INTERP_KERNEL::Exception(os.str().c_str());
235 if(_y_array->getNumberOfComponents()!=1)
237 std::ostringstream os; os << msg0 << 'Y' << msg2;
238 throw INTERP_KERNEL::Exception(os.str().c_str());
244 if(_z_array->getNbOfElems()<2)
246 std::ostringstream os; os << msg0 << 'Z' << msg1;
247 throw INTERP_KERNEL::Exception(os.str().c_str());
249 if(_z_array->getNumberOfComponents()!=1)
251 std::ostringstream os; os << msg0 << 'Z' << msg2;
252 throw INTERP_KERNEL::Exception(os.str().c_str());
257 void MEDCouplingCMesh::checkCoherency1(double eps) const throw(INTERP_KERNEL::Exception)
261 _x_array->checkMonotonic(true, eps);
263 _y_array->checkMonotonic(true, eps);
265 _z_array->checkMonotonic(true, eps);
268 void MEDCouplingCMesh::checkCoherency2(double eps) const throw(INTERP_KERNEL::Exception)
270 checkCoherency1(eps);
273 int MEDCouplingCMesh::getNumberOfCells() const
277 ret*=_x_array->getNbOfElems()-1;
279 ret*=_y_array->getNbOfElems()-1;
281 ret*=_z_array->getNbOfElems()-1;
285 int MEDCouplingCMesh::getNumberOfNodes() const
289 ret*=_x_array->getNbOfElems();
291 ret*=_y_array->getNbOfElems();
293 ret*=_z_array->getNbOfElems();
297 void MEDCouplingCMesh::getSplitCellValues(int *res) const
299 int spaceDim=getSpaceDimension();
300 for(int l=0;l<spaceDim;l++)
303 for(int p=0;p<spaceDim-l-1;p++)
304 val*=getCoordsAt(p)->getNbOfElems()-1;
305 res[spaceDim-l-1]=val;
309 void MEDCouplingCMesh::getSplitNodeValues(int *res) const
311 int spaceDim=getSpaceDimension();
312 for(int l=0;l<spaceDim;l++)
315 for(int p=0;p<spaceDim-l-1;p++)
316 val*=getCoordsAt(p)->getNbOfElems();
317 res[spaceDim-l-1]=val;
321 void MEDCouplingCMesh::getNodeGridStructure(int *res) const
323 int meshDim=getMeshDimension();
324 for(int i=0;i<meshDim;i++)
325 res[i]=getCoordsAt(i)->getNbOfElems();
328 int MEDCouplingCMesh::getSpaceDimension() const
340 int MEDCouplingCMesh::getMeshDimension() const
342 return getSpaceDimension();
345 void MEDCouplingCMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
348 int spaceDim=getSpaceDimension();
349 getSplitNodeValues(tmp);
350 const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
352 GetPosFromId(nodeId,spaceDim,tmp,tmp2);
353 for(int j=0;j<spaceDim;j++)
355 coo.push_back(tabs[j]->getConstPointer()[tmp2[j]]);
358 std::string MEDCouplingCMesh::simpleRepr() const
360 std::ostringstream ret;
361 ret << "Cartesian mesh with name : \"" << getName() << "\"\n";
362 ret << "Description of mesh : \"" << getDescription() << "\"\n";
364 double tt=getTime(tmpp1,tmpp2);
365 ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
366 ret << "Iteration : " << tmpp1 << " Order : " << tmpp2 << "\n";
367 ret << "Mesh and SpaceDimension dimension : " << getSpaceDimension() << "\n\nArrays :\n________\n\n";
370 ret << "X Array :\n";
371 _x_array->reprZipWithoutNameStream(ret);
375 ret << "Y Array :\n";
376 _y_array->reprZipWithoutNameStream(ret);
380 ret << "Z Array :\n";
381 _z_array->reprZipWithoutNameStream(ret);
386 std::string MEDCouplingCMesh::advancedRepr() const
391 const DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i) const throw(INTERP_KERNEL::Exception)
402 throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
406 DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i) throw(INTERP_KERNEL::Exception)
417 throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
421 void MEDCouplingCMesh::setCoordsAt(int i, const DataArrayDouble *arr) throw(INTERP_KERNEL::Exception)
424 arr->checkNbOfComps(1,"MEDCouplingCMesh::setCoordsAt");
425 DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
427 throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
428 if(arr!=*(thisArr[i]))
431 (*(thisArr[i]))->decrRef();
432 (*(thisArr[i]))=const_cast<DataArrayDouble *>(arr);
434 (*(thisArr[i]))->incrRef();
439 void MEDCouplingCMesh::setCoords(const DataArrayDouble *coordsX, const DataArrayDouble *coordsY, const DataArrayDouble *coordsZ)
442 coordsX->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsX");
444 coordsY->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsY");
446 coordsZ->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsZ");
449 _x_array=const_cast<DataArrayDouble *>(coordsX);
454 _y_array=const_cast<DataArrayDouble *>(coordsY);
459 _z_array=const_cast<DataArrayDouble *>(coordsZ);
465 void MEDCouplingCMesh::getBoundingBox(double *bbox) const
467 int dim=getSpaceDimension();
469 for (int idim=0; idim<dim; idim++)
471 const DataArrayDouble *c=getCoordsAt(idim);
474 const double *coords=c->getConstPointer();
475 int nb=c->getNbOfElems();
477 bbox[2*j+1]=coords[nb-1];
483 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
485 std::string name="MeasureOfMesh_";
487 int nbelem=getNumberOfCells();
488 MEDCouplingFieldDouble *field=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
489 field->setName(name.c_str());
490 DataArrayDouble* array=DataArrayDouble::New();
491 array->alloc(nbelem,1);
492 double *area_vol=array->getPointer();
493 field->setArray(array) ;
495 field->setMesh(const_cast<MEDCouplingCMesh *>(this));
496 field->synchronizeTimeWithMesh();
498 getSplitCellValues(tmp);
499 int dim=getSpaceDimension();
500 const double **thisArr=new const double *[dim];
501 const DataArrayDouble *thisArr2[3]={_x_array,_y_array,_z_array};
502 for(int i=0;i<dim;i++)
503 thisArr[i]=thisArr2[i]->getConstPointer();
504 for(int icell=0;icell<nbelem;icell++)
507 GetPosFromId(icell,dim,tmp,tmp2);
509 for(int i=0;i<dim;i++)
510 area_vol[icell]*=thisArr[i][tmp2[i]+1]-thisArr[i][tmp2[i]];
517 * not implemented yet !
519 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureFieldOnNode(bool isAbs) const
521 throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getMeasureFieldOnNode : not implemented yet !");
525 int MEDCouplingCMesh::getCellContainingPoint(const double *pos, double eps) const
527 int dim=getSpaceDimension();
530 for(int i=0;i<dim;i++)
532 const double *d=getCoordsAt(i)->getConstPointer();
533 int nbOfNodes=getCoordsAt(i)->getNbOfElems();
535 const double *w=std::find_if(d,d+nbOfNodes,std::bind2nd(std::greater_equal<double>(),ref));
536 int w2=(int)std::distance(d,w);
555 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
557 throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
560 void MEDCouplingCMesh::translate(const double *vector)
563 std::transform(_x_array->getConstPointer(),_x_array->getConstPointer()+_x_array->getNbOfElems(),
564 _x_array->getPointer(),std::bind2nd(std::plus<double>(),vector[0]));
566 std::transform(_y_array->getConstPointer(),_y_array->getConstPointer()+_y_array->getNbOfElems(),
567 _y_array->getPointer(),std::bind2nd(std::plus<double>(),vector[1]));
569 std::transform(_z_array->getConstPointer(),_z_array->getConstPointer()+_z_array->getNbOfElems(),
570 _z_array->getPointer(),std::bind2nd(std::plus<double>(),vector[2]));
573 void MEDCouplingCMesh::scale(const double *point, double factor)
577 DataArrayDouble *c=getCoordsAt(i);
580 double *coords=c->getPointer();
581 int lgth=c->getNbOfElems();
582 std::transform(coords,coords+lgth,coords,std::bind2nd(std::minus<double>(),point[i]));
583 std::transform(coords,coords+lgth,coords,std::bind2nd(std::multiplies<double>(),factor));
584 std::transform(coords,coords+lgth,coords,std::bind2nd(std::plus<double>(),point[i]));
591 MEDCouplingMesh *MEDCouplingCMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
593 //not implemented yet !
597 DataArrayDouble *MEDCouplingCMesh::getCoordinatesAndOwner() const
599 DataArrayDouble *ret=DataArrayDouble::New();
600 int spaceDim=getSpaceDimension();
601 int nbNodes=getNumberOfNodes();
602 ret->alloc(nbNodes,spaceDim);
603 double *pt=ret->getPointer();
605 getSplitNodeValues(tmp);
606 const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
607 const double *tabsPtr[3];
608 for(int j=0;j<spaceDim;j++)
610 tabsPtr[j]=tabs[j]->getConstPointer();
611 ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0).c_str());
614 for(int i=0;i<nbNodes;i++)
616 GetPosFromId(i,spaceDim,tmp,tmp2);
617 for(int j=0;j<spaceDim;j++)
618 pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
623 DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const
625 DataArrayDouble *ret=DataArrayDouble::New();
626 int spaceDim=getSpaceDimension();
627 int nbCells=getNumberOfCells();
628 ret->alloc(nbCells,spaceDim);
629 double *pt=ret->getPointer();
631 getSplitCellValues(tmp);
632 const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
633 std::vector<double> tabsPtr[3];
634 for(int j=0;j<spaceDim;j++)
636 int sz=tabs[j]->getNbOfElems()-1;
637 ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0).c_str());
638 const double *srcPtr=tabs[j]->getConstPointer();
639 tabsPtr[j].insert(tabsPtr[j].end(),srcPtr,srcPtr+sz);
640 std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),srcPtr+1,tabsPtr[j].begin(),std::plus<double>());
641 std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),tabsPtr[j].begin(),std::bind2nd(std::multiplies<double>(),0.5));
644 for(int i=0;i<nbCells;i++)
646 GetPosFromId(i,spaceDim,tmp,tmp2);
647 for(int j=0;j<spaceDim;j++)
648 pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
653 void MEDCouplingCMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
655 throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !");
658 void MEDCouplingCMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
661 double time=getTime(it,order);
664 littleStrings.clear();
665 littleStrings.push_back(getName());
666 littleStrings.push_back(getDescription());
667 littleStrings.push_back(getTimeUnit());
668 const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
675 val=thisArr[i]->getNumberOfTuples();
676 st=thisArr[i]->getInfoOnComponent(0);
678 tinyInfo.push_back(val);
679 littleStrings.push_back(st);
681 tinyInfo.push_back(it);
682 tinyInfo.push_back(order);
683 tinyInfoD.push_back(time);
686 void MEDCouplingCMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
696 void MEDCouplingCMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
698 a1=DataArrayInt::New();
700 const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
705 sz+=thisArr[i]->getNumberOfTuples();
707 a2=DataArrayDouble::New();
709 double *a2Ptr=a2->getPointer();
712 a2Ptr=std::copy(thisArr[i]->getConstPointer(),thisArr[i]->getConstPointer()+thisArr[i]->getNumberOfTuples(),a2Ptr);
715 void MEDCouplingCMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
716 const std::vector<std::string>& littleStrings)
718 setName(littleStrings[0].c_str());
719 setDescription(littleStrings[1].c_str());
720 setTimeUnit(littleStrings[2].c_str());
721 DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
722 const double *data=a2->getConstPointer();
727 (*(thisArr[i]))=DataArrayDouble::New();
728 (*(thisArr[i]))->alloc(tinyInfo[i],1);
729 (*(thisArr[i]))->setInfoOnComponent(0,littleStrings[i+3].c_str());
730 std::copy(data,data+tinyInfo[i],(*(thisArr[i]))->getPointer());
734 setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
737 void MEDCouplingCMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception)
739 std::ostringstream extent;
740 DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
744 { extent << "0 " << thisArr[i]->getNumberOfTuples()-1 << " "; }
746 { extent << "0 0 "; }
748 ofs << " <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\">\n";
749 ofs << " <Piece Extent=\"" << extent.str() << "\">\n";
750 ofs << " <PointData>\n" << pointData << std::endl;
751 ofs << " </PointData>\n";
752 ofs << " <CellData>\n" << cellData << std::endl;
753 ofs << " </CellData>\n";
754 ofs << " <Coordinates>\n";
758 thisArr[i]->writeVTK(ofs,8,"Array");
761 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=DataArrayDouble::New(); coo->alloc(1,1);
763 coo->writeVTK(ofs,8,"Array");
766 ofs << " </Coordinates>\n";
767 ofs << " </Piece>\n";
768 ofs << " </" << getVTKDataSetType() << ">\n";
771 std::string MEDCouplingCMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception)
773 return std::string("RectilinearGrid");