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