CDMATH/mesh/src/Field.cxx

   1 /*
   2  * field.cxx
   3  *
   4  *  Created on: 7 fevrier. 2012
   5  *      Author: CDMAT
   6  */
   7
   8 #include "Node.hxx"
   9 #include "Cell.hxx"
  10 #include "Field.hxx"
  11 #include "MEDFileMesh.hxx"
  12 #include "MEDLoader.hxx"
  13 #include "MEDCouplingUMesh.hxx"
  14 #include "MEDCouplingFieldDouble.hxx"
  15
  16 #include "CdmathException.hxx"
  17
  18 #include <fstream>
  19 #include <sstream>
  20 using namespace MEDCoupling;
  21 using namespace std;
  22
  23
  24 //----------------------------------------------------------------------
  25 Field::Field( EntityType typeField )
  26 //----------------------------------------------------------------------
  27 {
  28         _field=NULL;
  29         _typeField=typeField;
  30         _numberOfComponents=0;
  31 }
  32
  33 //----------------------------------------------------------------------
  34 Field::~Field( void )
  35 //----------------------------------------------------------------------
  36 {
  37         //std::cerr << "dtor Field, _field = " <<_field << std::endl;
  38         //if (_field) _field->decrRef();
  39 }
  40
  41
  42 Field::Field(const std::string fieldName, EntityType type, const Mesh& mesh, int numberOfComponents, double time)
  43 {
  44         _field = NULL;
  45         _mesh=Mesh(mesh);
  46         _typeField=type;
  47         _numberOfComponents=numberOfComponents;
  48         _time=time;
  49         _fieldName=fieldName;
  50
  51         buildFieldMemoryStructure();
  52 }
  53 void Field::buildFieldMemoryStructure()
  54 {
  55         MEDCouplingUMesh* mu=_mesh.getMEDCouplingMesh()->buildUnstructured();
  56         DataArrayDouble *array=DataArrayDouble::New();
  57         if (_typeField==CELLS)
  58         {
  59                 _field=MEDCouplingFieldDouble::New(ON_CELLS);
  60                 array->alloc(_mesh.getNumberOfCells(),_numberOfComponents);
  61                 _field->setMesh(mu);
  62         }else if(_typeField==NODES)
  63         {
  64                 _field=MEDCouplingFieldDouble::New(ON_NODES);
  65                 array->alloc(_mesh.getNumberOfNodes(),_numberOfComponents);
  66                 _field->setMesh(mu);
  67         }else if(_typeField==FACES)
  68         {
  69                 _field=MEDCouplingFieldDouble::New(ON_CELLS);
  70                 array->alloc(_mesh.getNumberOfFaces(),_numberOfComponents);
  71                 DataArrayInt *desc=DataArrayInt::New();
  72                 DataArrayInt *descI=DataArrayInt::New();
  73                 DataArrayInt *revDesc=DataArrayInt::New();
  74                 DataArrayInt *revDescI=DataArrayInt::New();
  75                 MEDCouplingUMesh *m3=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
  76                 _field->setMesh(m3);
  77                 desc->decrRef();
  78                 descI->decrRef();
  79                 revDesc->decrRef();
  80                 revDescI->decrRef();
  81                 m3->decrRef();
  82         }else
  83                 throw CdmathException("Type of Field::Field() is not compatible");
  84
  85         _field->setName(_fieldName.c_str()) ;
  86         _field->setArray(array);
  87         _field->setTime(_time,0,0);
  88         array->decrRef();
  89         mu->decrRef();
  90 }
  91
  92 Field::Field( const std::string filename, EntityType type,
  93                 const std::string & fieldName,
  94                 int iteration, int order, int meshLevel,
  95                 int numberOfComponents, double time)
  96 {
  97         _field = NULL;
  98         _mesh=Mesh(filename + ".med", meshLevel);
  99         _typeField=type;
 100         _numberOfComponents=numberOfComponents;
 101         _time=time;
 102         _fieldName=fieldName;
 103
 104         readFieldMed(filename, type, fieldName, iteration, order);
 105 }
 106
 107 Field::Field(const std::string meshFileName, EntityType type, const std::vector<double> Vconstant,
 108                 const std::string & fieldName, int meshLevel, double time )
 109 {
 110         _field = NULL;
 111         _mesh=Mesh(meshFileName + ".med", meshLevel);
 112         _typeField=type;
 113         _numberOfComponents=Vconstant.size();
 114         _time=time;
 115         _fieldName=fieldName;
 116
 117         buildFieldMemoryStructure();
 118
 119         int nbelem=_field->getNumberOfTuples();
 120         int nbcomp=_field->getNumberOfComponents() ;
 121
 122         for (int ielem=0 ; ielem<nbelem; ielem++)
 123                 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
 124                         _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
 125 }
 126 Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::string & fieldName, double time)
 127 {
 128         _field = NULL;
 129         _mesh=Mesh(M);
 130         _typeField=type;
 131         _numberOfComponents=Vconstant.size();
 132         _time=time;
 133         _fieldName=fieldName;
 134
 135         buildFieldMemoryStructure();
 136
 137         int nbelem=_field->getNumberOfTuples();
 138         int nbcomp=_field->getNumberOfComponents() ;
 139
 140         for (int ielem=0 ; ielem<nbelem; ielem++)
 141                 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
 142                         _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
 143 }
 144 Field::Field(const Mesh& M, EntityType type, const vector<double> Vconstant, const std::string & fieldName, double time)
 145 {
 146         _field = NULL;
 147         _mesh=Mesh(M);
 148         _typeField=type;
 149         _numberOfComponents=Vconstant.size();
 150         _time=time;
 151         _fieldName=fieldName;
 152
 153         buildFieldMemoryStructure();
 154
 155         int nbelem=_field->getNumberOfTuples();
 156         int nbcomp=_field->getNumberOfComponents() ;
 157
 158         for (int ielem=0 ; ielem<nbelem; ielem++)
 159                 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
 160                         _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
 161 }
 162 Field::Field( int nDim, const vector<double> Vconstant, EntityType type,
 163                 double xmin, double xmax, int nx, string leftSide, string rightSide,
 164                 double ymin, double ymax, int ny, string backSide, string frontSide,
 165                 double zmin, double zmax, int nz, string bottomSide, string topSide,
 166                 const std::string & fieldName, double time,double epsilon)
 167 {
 168         _field = NULL;
 169         _typeField=type;
 170         _numberOfComponents=Vconstant.size();
 171         _time=time;
 172         _fieldName=fieldName;
 173
 174         //Build mesh
 175         if(nDim==1){
 176                 _mesh=Mesh(xmin,xmax,nx);
 177         }
 178         else if(nDim==2)
 179                 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny);
 180         else if(nDim==3)
 181                 _mesh=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
 182         else{
 183                 cout<<"Field::Field: Space dimension nDim should be between 1 and 3"<<endl;
 184                 throw CdmathException("Space dimension nDim should be between 1 and 3");
 185         }
 186
 187         _mesh.setGroupAtPlan(xmax,0,epsilon,rightSide);
 188         _mesh.setGroupAtPlan(xmin,0,epsilon,leftSide);
 189         if(nDim>=2){
 190                 _mesh.setGroupAtPlan(ymax,1,epsilon,frontSide);
 191                 _mesh.setGroupAtPlan(ymin,1,epsilon,backSide);
 192         }
 193         if(nDim==3){
 194                 _mesh.setGroupAtPlan(zmax,2,epsilon,topSide);
 195                 _mesh.setGroupAtPlan(zmin,2,epsilon,bottomSide);
 196         }
 197
 198         // Build field
 199         buildFieldMemoryStructure();
 200
 201         int nbelem=_field->getNumberOfTuples();
 202         int nbcomp=_field->getNumberOfComponents() ;
 203
 204         for (int ielem=0 ; ielem<nbelem; ielem++)
 205                 for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
 206                         _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
 207 }
 208 Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos,
 209                 EntityType type, int direction, const std::string & fieldName, double time)
 210 {
 211         if  (VV_Right.getNumberOfRows()!=VV_Left.getNumberOfRows())
 212                 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
 213
 214         _field = NULL;
 215         _mesh=Mesh(M);
 216         _typeField=type;
 217         _numberOfComponents=VV_Left.getNumberOfRows();
 218         _time=time;
 219         _fieldName=fieldName;
 220
 221         // Build field
 222         buildFieldMemoryStructure();
 223
 224         int nbelem=_field->getNumberOfTuples();
 225         int nbcomp=_field->getNumberOfComponents() ;
 226         double component_value;
 227
 228         for (int j = 0; j < nbelem; j++) {
 229                 if(direction==0)
 230                         component_value=M.getCell(j).x();
 231                 else if(direction==1)
 232                         component_value=M.getCell(j).y();
 233                 else if(direction==2)
 234                         component_value=M.getCell(j).z();
 235                 else
 236                         throw CdmathException( "Field::Field: direction should be an integer between 0 and 2");
 237
 238                 for (int i=0; i< nbcomp; i++)
 239                         if (component_value< disc_pos )
 240                                 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Left[i];
 241                         else
 242                                 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Right[i];
 243         }
 244 }
 245 Field::Field( int nDim, const vector<double> VV_Left, vector<double> VV_Right,
 246                 double xstep, EntityType type,
 247                 double xmin, double xmax, int nx, string leftSide, string rightSide,
 248                 double ymin, double ymax, int ny, string backSide, string frontSide,
 249                 double zmin, double zmax, int nz, string bottomSide, string topSide,
 250                 int direction, const std::string & fieldName, double time, double epsilon)
 251 {
 252         if  (VV_Right.size()!=VV_Left.size())
 253                 throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
 254         Mesh M;
 255         if(nDim==1)
 256                 M=Mesh(xmin,xmax,nx);
 257         else if(nDim==2)
 258                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny);
 259         else if(nDim==3)
 260                 M=Mesh(xmin,xmax,nx,ymin,ymax,ny,zmin,zmax,nz);
 261         else
 262                 throw CdmathException("Field::Field : Space dimension nDim should be between 1 and 3");
 263
 264         M.setGroupAtPlan(xmax,0,epsilon,rightSide);
 265         M.setGroupAtPlan(xmin,0,epsilon,leftSide);
 266         if(nDim>=2){
 267                 M.setGroupAtPlan(ymax,1,epsilon,frontSide);
 268                 M.setGroupAtPlan(ymin,1,epsilon,backSide);
 269         }
 270         if(nDim==3){
 271                 M.setGroupAtPlan(zmax,2,epsilon,topSide);
 272                 M.setGroupAtPlan(zmin,2,epsilon,bottomSide);
 273         }
 274         Vector V_Left(VV_Left.size()), V_Right(VV_Right.size());
 275         for(int i=0;i<VV_Left.size(); i++){
 276                 V_Left(i)=VV_Left[i];
 277                 V_Right(i)=VV_Right[i];
 278         }
 279         Field(M, V_Left, V_Right, xstep,  type, direction, fieldName, time);
 280 }
 281
 282 Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius,
 283                 const Vector Center, EntityType type, const std::string & fieldName, double time)
 284 {
 285         if((Center.size()!=M.getSpaceDimension()) || (Vout.size() != Vin.size()) )
 286         {
 287                 cout<< "Vout.size()= "<<Vout.size() << ", Vin.size()= "<<Vin.size()<<", Center.size()="<<Center.size()<<", M.getSpaceDim= "<< M.getSpaceDimension()<<endl;
 288                 throw CdmathException("Field::Field : Vector size error");
 289         }
 290
 291         _field = NULL;
 292         _mesh=Mesh(M);
 293         _typeField=type;
 294         _numberOfComponents=Vout.size();
 295         _time=time;
 296         _fieldName=fieldName;
 297
 298         // Build field
 299         buildFieldMemoryStructure();
 300
 301         int nbelem=_field->getNumberOfTuples();
 302         int nbcomp=_field->getNumberOfComponents() ;
 303
 304         int spaceDim=M.getSpaceDimension();
 305         Vector currentPoint(spaceDim);
 306
 307         for(int i=0;i<nbelem;i++)
 308         {
 309                 currentPoint(0)=M.getCell(i).x();
 310                 if(spaceDim>1)
 311                 {
 312                         currentPoint(1)=M.getCell(i).y();
 313                         if(spaceDim>2)
 314                                 currentPoint(2)=M.getCell(i).z();
 315                 }
 316                 if((currentPoint-Center).norm()<radius)
 317                         for(int j=0;j<nbcomp;j++)
 318                                 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vin[j];
 319                 else
 320                         for(int j=0;j<nbcomp;j++)
 321                                 _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vout[j];
 322         }
 323 }
 324
 325 MEDCoupling::DataArrayDouble * Field::getArray(){
 326         return _field->getArray();
 327 }
 328
 329 void
 330 Field::readFieldMed( const std::string & fileNameRadical,
 331                 EntityType type,
 332                 const std::string & fieldName,
 333                 int iteration,
 334                 int order)
 335 {
 336         /**
 337          * Reads the file fileNameRadical.med and creates a Field from it.
 338          */
 339         std::string completeFileName = fileNameRadical + ".med";
 340         std::vector<std::string> fieldNames = MEDCoupling::GetAllFieldNames(completeFileName);
 341         size_t iField = 0;
 342         std::string attributedFieldName;
 343         _field = NULL;
 344
 345         // Get the name of the right field that we will attribute to the Field.
 346         if (fieldName == "") {
 347                 if (fieldNames.size() > 0)
 348                         attributedFieldName = fieldNames[0];
 349                 else {
 350                         std::ostringstream message;
 351                         message << "No field in file " << completeFileName;
 352                         throw CdmathException(message.str().c_str());
 353                 }
 354         }
 355         else {
 356                 for (; iField < fieldNames.size(); iField++)
 357                         if (fieldName == fieldNames[iField]) break;
 358
 359                 if (iField < fieldNames.size())
 360                         attributedFieldName = fieldName;
 361                 else {
 362                         std::ostringstream message;
 363                         message << "No field named " << fieldName << " in file " << completeFileName;
 364                         throw CdmathException(message.str().c_str());
 365                 }
 366         }
 367
 368         // Get the name of the right mesh that we will attribute to the Field.
 369         std::vector<std::string> meshNames
 370         = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName);
 371         if (meshNames.size() == 0) {
 372                 std::ostringstream message;
 373                 message << "No mesh associated to " << fieldName
 374                                 << " in file " << completeFileName;
 375                 throw CdmathException(message.str().c_str());
 376         }
 377         std::string attributedMeshName = meshNames[0];
 378
 379         // Create Field.
 380         MEDCoupling::TypeOfField medFieldType[3] = { ON_CELLS, ON_NODES, ON_CELLS };
 381         switch (type) {
 382         case CELLS:
 383                 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
 384                                         MEDCoupling::ReadFieldCell( completeFileName,
 385                                         attributedMeshName, 0,
 386                                         attributedFieldName, iteration, order) );
 387                 break;
 388         case NODES:
 389                 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
 390                                         MEDCoupling::ReadFieldNode( completeFileName,
 391                                         attributedMeshName, 0,
 392                                         attributedFieldName, iteration, order) );
 393                 break;
 394         case FACES:
 395                 _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
 396                                         MEDCoupling::ReadFieldCell( completeFileName,
 397                                         attributedMeshName, -1,
 398                                         attributedFieldName, iteration, order) );
 399                 break;
 400         }
 401 }
 402
 403
 404 DoubleTab
 405 Field::getNormEuclidean() const
 406 {
 407         DoubleTab norm(getNumberOfElements(),_field->magnitude()->getArray()->getConstPointer());
 408         return norm;
 409 }
 410
 411 double
 412 Field::max() const
 413 {
 414         if( getNumberOfComponents() !=1)
 415                 throw CdmathException("double Field::max() : field should have a single component in order to extract maximum value");
 416
 417         double result=0;
 418         for(int i=0; i<getNumberOfElements() ; i++)
 419                 if( result < (*this)(i,0))
 420                         result = (*this)(i,0);
 421
 422         return result;
 423 }
 424
 425 double
 426 Field::min() const
 427 {
 428         if( getNumberOfComponents() !=1)
 429                 throw CdmathException("double Field::min() : field should have a single component in order to extract minimum value");
 430
 431         double result=1e100;
 432         for(int i=0; i<getNumberOfElements() ; i++)
 433                 if( result > (*this)(i,0))
 434                         result = (*this)(i,0);
 435
 436         return result;
 437 }
 438
 439 string
 440 Field::getInfoOnComponent(int icomp) const
 441 {
 442         return _field->getArray()->getInfoOnComponent(icomp);
 443 }
 444
 445 void
 446 Field::setInfoOnComponent(int icomp, string nameCompo)
 447 {
 448         _field.retn()->getArray()->setInfoOnComponent(icomp,nameCompo);
 449 }
 450
 451 //----------------------------------------------------------------------
 452 Field::Field( const Field & f )
 453 //----------------------------------------------------------------------
 454 {
 455         _mesh=f.getMesh() ;
 456         MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
 457         _field=f1;
 458         _typeField=f.getTypeOfField();
 459 }
 460
 461 //----------------------------------------------------------------------
 462 MCAuto<MEDCouplingFieldDouble>
 463 Field::getField ( void )  const
 464 //----------------------------------------------------------------------
 465 {
 466         return _field ;
 467 }
 468
 469 //----------------------------------------------------------------------
 470 void
 471 Field::setFieldByMEDCouplingFieldDouble ( const MEDCouplingFieldDouble* field )
 472 //----------------------------------------------------------------------
 473 {
 474         MCAuto<MEDCouplingFieldDouble> ff=field->deepCopy();
 475         _field=ff;
 476 }
 477
 478 //----------------------------------------------------------------------
 479 void
 480 Field::setFieldByDataArrayDouble ( const DataArrayDouble* array )
 481 //----------------------------------------------------------------------
 482 {
 483         _field->setArray(const_cast<DataArrayDouble*>(array));
 484 }
 485
 486 //----------------------------------------------------------------------
 487 Vector
 488 Field::integral ( ) const
 489 //----------------------------------------------------------------------
 490 {
 491         int nbComp=_field->getNumberOfComponents();
 492         double res[nbComp];//Pointer containing the inegral of each component
 493         _field->integral(false,res);
 494         Vector result(nbComp);//Vector containing the inegral of each component
 495
 496         for(int i=0; i<nbComp ; i++)
 497                 result(i)=res[i];
 498
 499         return result;
 500 }
 501
 502 //----------------------------------------------------------------------
 503 double
 504 Field::integral ( int compId ) const
 505 //----------------------------------------------------------------------
 506 {
 507         return _field->integral(compId, false);
 508 }
 509
 510 //----------------------------------------------------------------------
 511 Vector
 512 Field::normL1 ( ) const
 513 //----------------------------------------------------------------------
 514 {
 515         int nbComp=_field->getNumberOfComponents();
 516         double res[nbComp];//Pointer containing the L1 norm of each component
 517         _field->normL1(res);
 518         Vector result(nbComp);//Vector containing the L1 norm of each component
 519
 520         for(int i=0; i<nbComp ; i++)
 521                 result(i)=res[i];
 522
 523         return result;
 524 }
 525
 526 //----------------------------------------------------------------------
 527 Vector
 528 Field::normL2 ( ) const
 529 //----------------------------------------------------------------------
 530 {
 531         int nbComp=_field->getNumberOfComponents();
 532         double res[nbComp];//Pointer containing the L2 norm of each component
 533         _field->normL2(res);
 534         Vector result(nbComp);//Vector containing the L2 norm of each component
 535
 536         for(int i=0; i<nbComp ; i++)
 537                 result(i)=res[i];
 538
 539         return result;
 540 }
 541
 542 //----------------------------------------------------------------------
 543 Vector
 544 Field::normMax ( ) const
 545 //----------------------------------------------------------------------
 546 {
 547         int nbComp=_field->getNumberOfComponents();
 548         int nbElems=getNumberOfElements();
 549
 550         Vector result(nbComp);//Vector containing the Linfinity norm of each component
 551
 552         for(int i=0; i<nbElems ; i++)
 553         for(int j=0; j<nbComp ; j++)
 554             if(fabs((*this)(i,j))>result(j))
 555                 result(j)=fabs((*this)(i,j));
 556
 557         return result;
 558 }
 559
 560 //----------------------------------------------------------------------
 561 double&
 562 Field::operator() ( int ielem )
 563 //----------------------------------------------------------------------
 564 {
 565         if(ielem>_field->getNumberOfTuples() || ielem<0)
 566                 throw CdmathException("double& Field::operator(ielem) : ielem>number of values !");
 567         return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
 568 }
 569
 570 //----------------------------------------------------------------------
 571 double&
 572 Field::operator[] ( int ielem )
 573 //----------------------------------------------------------------------
 574 {
 575         if(ielem>_field->getNumberOfTuples() || ielem<0)
 576                 throw CdmathException("double& Field::operator[ielem] : ielem>number of values !");
 577         return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
 578 }
 579
 580 //----------------------------------------------------------------------
 581 double
 582 Field::operator() ( int ielem ) const
 583 //----------------------------------------------------------------------
 584 {
 585         if(ielem>_field->getNumberOfTuples() || ielem<0)
 586                 throw CdmathException("double Field::operator(ielem) : ielem>number of values !");
 587         return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
 588 }
 589
 590 //----------------------------------------------------------------------
 591 double
 592 Field::operator[] ( int ielem ) const
 593 //----------------------------------------------------------------------
 594 {
 595         if(ielem>_field->getNumberOfTuples() || ielem<0)
 596                 throw CdmathException("double Field::operator[ielem] : ielem>number of values !");
 597         return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
 598 }
 599
 600 //----------------------------------------------------------------------
 601 double&
 602 Field::operator() ( int ielem, int jcomp )
 603 //----------------------------------------------------------------------
 604 {
 605         if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
 606                 throw CdmathException("double& Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
 607         return _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()];
 608 }
 609
 610 //----------------------------------------------------------------------
 611 double
 612 Field::operator() (  int ielem, int jcomp ) const
 613 //----------------------------------------------------------------------
 614 {
 615         if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
 616                 throw CdmathException("double Field::operator(  int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
 617         return _field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()];
 618 }
 619
 620 //----------------------------------------------------------------------
 621 void
 622 Field::setTime ( double time, int iter )
 623 //----------------------------------------------------------------------
 624 {
 625         _field->setTime(time,iter,0.0);
 626 }
 627 //----------------------------------------------------------------------
 628 double
 629 Field::getTime ( void ) const
 630 //----------------------------------------------------------------------
 631 {
 632         int a,b;
 633         return _field->getTime(a,b);
 634 }
 635
 636 //----------------------------------------------------------------------
 637 int
 638 Field::getNumberOfElements ( void ) const
 639 //----------------------------------------------------------------------
 640 {
 641         return _field->getNumberOfTuples() ;
 642 }
 643
 644 int
 645 Field::getSpaceDimension( void ) const
 646 {
 647         return _mesh.getSpaceDimension() ;
 648 }
 649
 650 //----------------------------------------------------------------------
 651 int
 652 Field::getNumberOfComponents ( void ) const
 653 //----------------------------------------------------------------------
 654 {
 655         return _field->getNumberOfComponents() ;
 656 }
 657
 658 //----------------------------------------------------------------------
 659 const double*
 660 Field::getValues ( void ) const
 661 //----------------------------------------------------------------------
 662 {
 663         return _field->getArray()->getConstPointer() ;
 664 }
 665
 666 //----------------------------------------------------------------------
 667 const string
 668 Field::getName ( void ) const
 669 //----------------------------------------------------------------------
 670 {
 671         return _field->getName() ;
 672 }
 673
 674 //----------------------------------------------------------------------
 675 const Mesh&
 676 Field::getMesh ( void ) const
 677 //----------------------------------------------------------------------
 678 {
 679         return _mesh ;
 680 }
 681
 682 //----------------------------------------------------------------------
 683 EntityType
 684 Field::getTypeOfField ( void ) const
 685 //----------------------------------------------------------------------
 686 {
 687         return _typeField;
 688 }
 689
 690 //----------------------------------------------------------------------
 691 void
 692 Field::setName ( const string fieldName )
 693 //----------------------------------------------------------------------
 694 {
 695         _field->setName(fieldName.c_str()) ;
 696 }
 697
 698
 699 //----------------------------------------------------------------------
 700 Field
 701 Field::operator+ ( const Field& f ) const
 702 //----------------------------------------------------------------------
 703 {
 704         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 705         //throw CdmathException("Field::operator+ : Field addition requires identical meshes");
 706         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 707         if(f.getTypeOfField() != getTypeOfField())
 708                 throw CdmathException("Field::operator+ : Field addition requires identical field types (CELLS, NODES or FACES");
 709         if(f.getNumberOfComponents() != getNumberOfComponents())
 710                 throw CdmathException("Field::operator+ : Field addition requires identical number of components");
 711
 712         Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
 713         int nbComp=f.getNumberOfComponents();
 714         int nbElem=f.getNumberOfElements();
 715         for (int ielem=0 ; ielem<nbElem; ielem++)
 716                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
 717                         fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp);
 718         return fres;
 719 }
 720
 721 //----------------------------------------------------------------------
 722 Field
 723 Field::operator- ( const Field& f ) const
 724 //----------------------------------------------------------------------
 725 {
 726         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 727         //throw CdmathException("Field::operator- : Field subtraction requires identical meshes");
 728         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 729         if(f.getTypeOfField() != getTypeOfField())
 730                 throw CdmathException("Field::operator- : Field subtraction requires identical field types (CELLS, NODES or FACES");
 731         if(f.getNumberOfComponents() != getNumberOfComponents())
 732                 throw CdmathException("Field::operator- : Field subtraction requires identical number of components");
 733
 734         Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
 735         int nbComp=f.getNumberOfComponents();
 736         int nbElem=f.getNumberOfElements();
 737         for (int ielem=0 ; ielem<nbElem; ielem++)
 738                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
 739                         fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp);
 740         return fres;
 741 }
 742
 743 //----------------------------------------------------------------------
 744 const Field&
 745 Field::operator= ( const Field& f )
 746 //----------------------------------------------------------------------
 747 {
 748         _mesh=f.getMesh() ;
 749         _typeField=f.getTypeOfField() ;
 750         _numberOfComponents=f.getNumberOfComponents();
 751         _time=f.getTime();
 752         _fieldName=f.getName();
 753         MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
 754         _field=f1;
 755         return *this;
 756 }
 757
 758 Field Field::deepCopy( ) const
 759 {
 760     Field F(getName(), getTypeOfField(), getMesh(), getNumberOfComponents(), getTime()) ;
 761         MCAuto<MEDCouplingFieldDouble> f1=getField()->deepCopy();
 762         F.setFieldByMEDCouplingFieldDouble(f1);
 763
 764     return F;
 765 }
 766
 767
 768 //----------------------------------------------------------------------
 769 const Field&
 770 Field::operator+= ( const Field& f )
 771 //----------------------------------------------------------------------
 772 {
 773         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 774         //throw CdmathException("Field::operator+= : Field addition requires identical meshes");
 775         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 776         if(f.getTypeOfField() != getTypeOfField())
 777                 throw CdmathException("Field::operator+= : Field addition requires identical field types (CELLS, NODES or FACES");
 778         if(f.getNumberOfComponents() != getNumberOfComponents())
 779                 throw CdmathException("Field::operator+= : Field addition requires identical number of components");
 780
 781         _field->setMesh(f.getField()->getMesh());
 782         (*_field)+=(*f.getField());
 783         return *this;
 784 }
 785
 786 //----------------------------------------------------------------------
 787 const Field&
 788 Field::operator-= ( const Field& f )
 789 //----------------------------------------------------------------------
 790 {
 791         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 792         //throw CdmathException("Field::operator-= : Field subtraction requires identical meshes");
 793         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 794         if(f.getTypeOfField() != getTypeOfField())
 795                 throw CdmathException("Field::operator-= : Field subtraction requires identical field types (CELLS, NODES or FACES");
 796         if(f.getNumberOfComponents() != getNumberOfComponents())
 797                 throw CdmathException("Field::operator-= : Field subtraction requires identical number of components");
 798
 799         _field->setMesh(f.getField()->getMesh());
 800         (*_field)-=(*f.getField());
 801         return *this;
 802 }
 803
 804 //----------------------------------------------------------------------
 805 const Field&
 806 Field::operator*= ( double s )
 807 //----------------------------------------------------------------------
 808 {
 809         int nbComp=getNumberOfComponents();
 810         int nbElem=getNumberOfElements();
 811         for (int i=0 ; i<nbComp ; i++)
 812                 for (int j=0 ; j<nbElem; j++)
 813                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s;
 814         return *this;
 815 }
 816
 817 //----------------------------------------------------------------------
 818 const Field&
 819 Field::operator/= ( double s )
 820 //----------------------------------------------------------------------
 821 {
 822         int nbComp=getNumberOfComponents();
 823         int nbElem=getNumberOfElements();
 824         for (int i=0 ; i<nbComp ; i++)
 825                 for (int j=0 ; j<nbElem; j++)
 826                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s;
 827         return *this;
 828 }
 829
 830 //----------------------------------------------------------------------
 831 const Field&
 832 Field::operator-= ( double s )
 833 //----------------------------------------------------------------------
 834 {
 835         int nbComp=getNumberOfComponents();
 836         int nbElem=getNumberOfElements();
 837         for (int i=0 ; i<nbComp ; i++)
 838                 for (int j=0 ; j<nbElem; j++)
 839                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s;
 840         return *this;
 841 }
 842
 843 //----------------------------------------------------------------------
 844 const Field&
 845 Field::operator+= ( double s )
 846 //----------------------------------------------------------------------
 847 {
 848         int nbComp=getNumberOfComponents();
 849         int nbElem=getNumberOfElements();
 850         for (int i=0 ; i<nbComp ; i++)
 851                 for (int j=0 ; j<nbElem; j++)
 852                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s;
 853         return *this;
 854 }
 855
 856 //----------------------------------------------------------------------
 857 void
 858 Field::writeVTK (std::string fileName, bool fromScratch) const
 859 //----------------------------------------------------------------------
 860 {
 861         string fname=fileName+".pvd";
 862         int iter,order;
 863         double time=_field->getTime(iter,order);
 864
 865         if (fromScratch)
 866         {
 867                 ofstream file(fname.c_str()) ;
 868                 file << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\"><Collection>\n" ;
 869                 ostringstream numfile;
 870                 numfile << iter ;
 871                 string filetmp=fileName+"_";
 872                 filetmp=filetmp+numfile.str();
 873                 string ret=_field->writeVTK(filetmp.c_str()) ;
 874                 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
 875                 file << "</Collection></VTKFile>\n" ;
 876                 file.close() ;
 877         }
 878         else
 879         {
 880                 ifstream file1(fname.c_str()) ;
 881                 string contenus;
 882                 getline(file1, contenus, '\0');
 883                 string to_remove="</Collection></VTKFile>";
 884                 size_t m = contenus.find(to_remove);
 885                 size_t n = contenus.find_first_of("\n", m + to_remove.length());
 886                 contenus.erase(m, n - m + 1);
 887                 file1.close() ;
 888                 ofstream file(fname.c_str()) ;
 889                 file << contenus ;
 890                 ostringstream numfile;
 891                 numfile << iter ;
 892                 string filetmp=fileName+"_";
 893                 filetmp=filetmp+numfile.str();
 894                 string ret=_field->writeVTK(filetmp.c_str()) ;
 895                 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
 896                 file << "</Collection></VTKFile>\n" ;
 897                 file.close() ;
 898         }
 899 }
 900
 901 //----------------------------------------------------------------------
 902 void
 903 Field::writeCSV ( const std::string fileName ) const
 904 //----------------------------------------------------------------------
 905 {
 906         int iter,order;
 907         double time=_field->getTime(iter,order);
 908
 909         ostringstream numfile;
 910         numfile << iter ;
 911         string filetmp=fileName+"_";
 912         filetmp=filetmp+numfile.str();
 913         filetmp=filetmp+".csv";
 914         ofstream file(filetmp.c_str()) ;
 915         int dim=_mesh.getSpaceDimension();
 916         int nbElements;
 917         if (getTypeOfField()==CELLS)
 918                 nbElements=_mesh.getNumberOfCells();
 919         else
 920                 nbElements=_mesh.getNumberOfNodes();
 921
 922         if (dim==1)
 923         {
 924                 int nbCompo=getNumberOfComponents();
 925                 if (nbCompo==1)
 926                         file << "x " << _field->getName() << endl;
 927                 else if (nbCompo>1)
 928                 {
 929                         file << "x";
 930                         for (int i=0;i<nbCompo;i++)
 931                                 file << " " << _field->getName() << " Compo " << i+1 << " "<< getInfoOnComponent(i);
 932                         file << endl;
 933                 }
 934                 for (int i=0;i<nbElements;i++)
 935                 {
 936                         if (getTypeOfField()==CELLS)
 937                                 file << _mesh.getCell(i).x() ;
 938                         else
 939                                 file << _mesh.getNode(i).x() ;
 940                         for (int j=0;j<nbCompo;j++)
 941                                 file << " " << getValues()[j+i*nbCompo] ;
 942                         file << endl;
 943                 }
 944         }else if (dim==2)
 945         {
 946                 int nbCompo=getNumberOfComponents();
 947                 if (nbCompo==1)
 948                         file << "x y " << _field->getName() << endl;
 949                 else if (nbCompo>1)
 950                 {
 951                         file << "x y";
 952                         for (int i=0;i<nbCompo;i++)
 953                                 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
 954                         file << endl;
 955                 }
 956                 for (int i=0;i<nbElements;i++)
 957                 {
 958                         if (getTypeOfField()==CELLS)
 959                                 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() ;
 960                         else
 961                                 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() ;
 962                         for (int j=0;j<nbCompo;j++)
 963                                 file << " " << getValues()[j+i*nbCompo] ;
 964                         file << endl;
 965                 }
 966         }else
 967         {
 968                 int nbCompo=getNumberOfComponents();
 969                 if (nbCompo==1)
 970                         file << "x y z " << _field->getName() << endl;
 971                 else if (nbCompo>1)
 972                 {
 973                         file << "x y z";
 974                         for (int i=0;i<nbCompo;i++)
 975                                 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
 976                         file << endl;
 977                 }
 978                 for (int i=0;i<nbElements;i++)
 979                 {
 980                         if (getTypeOfField()==CELLS)
 981                                 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() << " " << _mesh.getCell(i).z();
 982                         else
 983                                 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() << " " << _mesh.getNode(i).z();
 984                         for (int j=0;j<nbCompo;j++)
 985                                 file << " " << getValues()[j+i*nbCompo] ;
 986                         file << endl;
 987                 }
 988         }
 989         file.close() ;
 990 }
 991
 992 //----------------------------------------------------------------------
 993 void
 994 Field::writeMED ( const std::string fileName, bool fromScratch) const
 995 //----------------------------------------------------------------------
 996 {
 997         string fname=fileName+".med";
 998         if (fromScratch)
 999                 MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
1000         else
1001                 MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
1002 }
1003
1004 Field
1005 operator* (double value , const Field& field )
1006 {
1007         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1008         int nbComp=field.getNumberOfComponents();
1009         int nbElem=field.getNumberOfElements();
1010         for (int ielem=0 ; ielem<nbElem; ielem++)
1011                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1012                         fres(ielem, jcomp)=value*field(ielem, jcomp);
1013         return fres;
1014 }
1015
1016 Field
1017 operator* (const Field& field, double value )
1018 {
1019         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1020         int nbComp=field.getNumberOfComponents();
1021         int nbElem=field.getNumberOfElements();
1022         for (int ielem=0 ; ielem<nbElem; ielem++)
1023                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1024                         fres(ielem, jcomp)=value*field(ielem, jcomp);
1025         return fres;
1026 }
1027
1028 Field operator/ (const Field& field, double value)
1029                                 {
1030         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1031         int nbComp=field.getNumberOfComponents();
1032         int nbElem=field.getNumberOfElements();
1033         for (int ielem=0 ; ielem<nbElem; ielem++)
1034                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1035                         fres(ielem, jcomp)=field(ielem, jcomp)/value;
1036         return fres;
1037                                 }
1038
1039 Vector
1040 Field::getValuesOnAllComponents(int elem) const
1041 {
1042         Vector v(getNumberOfComponents());
1043         for(int i=0;i<getNumberOfComponents();i++)
1044                 v[i]=(*this)(elem,i);
1045         return v;
1046 }
1047
1048 Vector
1049 Field::getValuesOnComponent(int compo) const
1050 {
1051         Vector v(getNumberOfElements());
1052         for(int i=0;i<getNumberOfElements();i++)
1053                 v[i]=(*this)(i,compo);
1054         return v;
1055 }
1056
1057 std::ostream& operator<<(std::ostream& out, const Field& field )
1058 {
1059         cout << "Field " << field.getName() << " : " << endl ;
1060         out<< field.getField().retn()->getArray()->repr();
1061         return out;
1062 }