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(int component) const
 413 {
 414         if( component >= getNumberOfComponents() )
 415                 throw CdmathException("double Field::max() : component number should be smaller than field number of components");
 416
 417         double result=-1e100;
 418         for(int i=0; i<getNumberOfElements() ; i++)
 419                 if( result < (*this)(i,component))
 420                         result = (*this)(i,component);
 421
 422         return result;
 423 }
 424
 425 double
 426 Field::min(int component) const
 427 {
 428         if( component >= getNumberOfComponents() )
 429                 throw CdmathException("double Field::min() : component number should be smaller than field number of components");
 430
 431         double result=1e100;
 432         for(int i=0; i<getNumberOfElements() ; i++)
 433                 if( result > (*this)(i,component))
 434                         result = (*this)(i,component);
 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 Vector
 561 Field::componentMax(Vector & Indices) const
 562 {
 563         int nbComp=_field->getNumberOfComponents();
 564         int nbElems=getNumberOfElements();
 565
 566         Vector result(nbComp);//Vector containing the Linfinity norm of each component
 567
 568         for(int i=0; i<nbElems ; i++)
 569         for(int j=0; j<nbComp ; j++)
 570             if(fabs((*this)(i,j))>result(j))
 571             {
 572                 result(j)=fabs((*this)(i,j));
 573                 Indices(j)=i;
 574             }
 575         return result;
 576 }
 577
 578 //----------------------------------------------------------------------
 579 double&
 580 Field::operator() ( int ielem )
 581 //----------------------------------------------------------------------
 582 {
 583         if(ielem>_field->getNumberOfTuples() || ielem<0)
 584                 throw CdmathException("double& Field::operator(ielem) : ielem>number of values !");
 585         return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
 586 }
 587
 588 //----------------------------------------------------------------------
 589 double&
 590 Field::operator[] ( int ielem )
 591 //----------------------------------------------------------------------
 592 {
 593         if(ielem>_field->getNumberOfTuples() || ielem<0)
 594                 throw CdmathException("double& Field::operator[ielem] : ielem>number of values !");
 595         return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
 596 }
 597
 598 //----------------------------------------------------------------------
 599 double
 600 Field::operator() ( int ielem ) const
 601 //----------------------------------------------------------------------
 602 {
 603         if(ielem>_field->getNumberOfTuples() || ielem<0)
 604                 throw CdmathException("double Field::operator(ielem) : ielem>number of values !");
 605         return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
 606 }
 607
 608 //----------------------------------------------------------------------
 609 double
 610 Field::operator[] ( int ielem ) const
 611 //----------------------------------------------------------------------
 612 {
 613         if(ielem>_field->getNumberOfTuples() || ielem<0)
 614                 throw CdmathException("double Field::operator[ielem] : ielem>number of values !");
 615         return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
 616 }
 617
 618 //----------------------------------------------------------------------
 619 double&
 620 Field::operator() ( int ielem, int jcomp )
 621 //----------------------------------------------------------------------
 622 {
 623         if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
 624                 throw CdmathException("double& Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
 625         return _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()];
 626 }
 627
 628 //----------------------------------------------------------------------
 629 double
 630 Field::operator() (  int ielem, int jcomp ) const
 631 //----------------------------------------------------------------------
 632 {
 633         if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
 634                 throw CdmathException("double Field::operator(  int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
 635         return _field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()];
 636 }
 637
 638 //----------------------------------------------------------------------
 639 void
 640 Field::setTime ( double time, int iter )
 641 //----------------------------------------------------------------------
 642 {
 643         _field->setTime(time,iter,0.0);
 644 }
 645 //----------------------------------------------------------------------
 646 double
 647 Field::getTime ( void ) const
 648 //----------------------------------------------------------------------
 649 {
 650         int a,b;
 651         return _field->getTime(a,b);
 652 }
 653
 654 //----------------------------------------------------------------------
 655 int
 656 Field::getNumberOfElements ( void ) const
 657 //----------------------------------------------------------------------
 658 {
 659         return _field->getNumberOfTuples() ;
 660 }
 661
 662 int
 663 Field::getSpaceDimension( void ) const
 664 {
 665         return _mesh.getSpaceDimension() ;
 666 }
 667
 668 //----------------------------------------------------------------------
 669 int
 670 Field::getNumberOfComponents ( void ) const
 671 //----------------------------------------------------------------------
 672 {
 673         return _field->getNumberOfComponents() ;
 674 }
 675
 676 //----------------------------------------------------------------------
 677 const double*
 678 Field::getValues ( void ) const
 679 //----------------------------------------------------------------------
 680 {
 681         return _field->getArray()->getConstPointer() ;
 682 }
 683
 684 //----------------------------------------------------------------------
 685 const string
 686 Field::getName ( void ) const
 687 //----------------------------------------------------------------------
 688 {
 689         return _field->getName() ;
 690 }
 691
 692 //----------------------------------------------------------------------
 693 const Mesh&
 694 Field::getMesh ( void ) const
 695 //----------------------------------------------------------------------
 696 {
 697         return _mesh ;
 698 }
 699
 700 //----------------------------------------------------------------------
 701 EntityType
 702 Field::getTypeOfField ( void ) const
 703 //----------------------------------------------------------------------
 704 {
 705         return _typeField;
 706 }
 707
 708 //----------------------------------------------------------------------
 709 void
 710 Field::setName ( const string fieldName )
 711 //----------------------------------------------------------------------
 712 {
 713         _field->setName(fieldName.c_str()) ;
 714 }
 715
 716
 717 //----------------------------------------------------------------------
 718 Field
 719 Field::operator+ ( const Field& f ) const
 720 //----------------------------------------------------------------------
 721 {
 722         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 723         //throw CdmathException("Field::operator+ : Field addition requires identical meshes");
 724         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 725         if(f.getTypeOfField() != getTypeOfField())
 726                 throw CdmathException("Field::operator+ : Field addition requires identical field types (CELLS, NODES or FACES");
 727         if(f.getNumberOfComponents() != getNumberOfComponents())
 728                 throw CdmathException("Field::operator+ : Field addition requires identical number of components");
 729
 730         Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
 731         int nbComp=f.getNumberOfComponents();
 732         int nbElem=f.getNumberOfElements();
 733         for (int ielem=0 ; ielem<nbElem; ielem++)
 734                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
 735                         fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp);
 736         return fres;
 737 }
 738
 739 //----------------------------------------------------------------------
 740 Field
 741 Field::operator- ( const Field& f ) const
 742 //----------------------------------------------------------------------
 743 {
 744         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 745         //throw CdmathException("Field::operator- : Field subtraction requires identical meshes");
 746         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 747         if(f.getTypeOfField() != getTypeOfField())
 748                 throw CdmathException("Field::operator- : Field subtraction requires identical field types (CELLS, NODES or FACES");
 749         if(f.getNumberOfComponents() != getNumberOfComponents())
 750                 throw CdmathException("Field::operator- : Field subtraction requires identical number of components");
 751
 752         Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
 753         int nbComp=f.getNumberOfComponents();
 754         int nbElem=f.getNumberOfElements();
 755         for (int ielem=0 ; ielem<nbElem; ielem++)
 756                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
 757                         fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp);
 758         return fres;
 759 }
 760
 761 //----------------------------------------------------------------------
 762 const Field&
 763 Field::operator= ( const Field& f )
 764 //----------------------------------------------------------------------
 765 {
 766         _mesh=f.getMesh() ;
 767         _typeField=f.getTypeOfField() ;
 768         _numberOfComponents=f.getNumberOfComponents();
 769         _time=f.getTime();
 770         _fieldName=f.getName();
 771         MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
 772         _field=f1;
 773         return *this;
 774 }
 775
 776 Field Field::deepCopy( ) const
 777 {
 778     Field F(getName(), getTypeOfField(), getMesh(), getNumberOfComponents(), getTime()) ;
 779         MCAuto<MEDCouplingFieldDouble> f1=getField()->deepCopy();
 780         F.setFieldByMEDCouplingFieldDouble(f1);
 781
 782     return F;
 783 }
 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 addition requires identical meshes");
 793         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 794         if(f.getTypeOfField() != getTypeOfField())
 795                 throw CdmathException("Field::operator+= : Field addition requires identical field types (CELLS, NODES or FACES");
 796         if(f.getNumberOfComponents() != getNumberOfComponents())
 797                 throw CdmathException("Field::operator+= : Field addition 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-= ( const Field& f )
 807 //----------------------------------------------------------------------
 808 {
 809         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 810         //throw CdmathException("Field::operator-= : Field subtraction requires identical meshes");
 811         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 812         if(f.getTypeOfField() != getTypeOfField())
 813                 throw CdmathException("Field::operator-= : Field subtraction requires identical field types (CELLS, NODES or FACES");
 814         if(f.getNumberOfComponents() != getNumberOfComponents())
 815                 throw CdmathException("Field::operator-= : Field subtraction requires identical number of components");
 816
 817         _field->setMesh(f.getField()->getMesh());
 818         (*_field)-=(*f.getField());
 819         return *this;
 820 }
 821
 822 //----------------------------------------------------------------------
 823 const Field&
 824 Field::operator*= ( double s )
 825 //----------------------------------------------------------------------
 826 {
 827         int nbComp=getNumberOfComponents();
 828         int nbElem=getNumberOfElements();
 829         for (int i=0 ; i<nbComp ; i++)
 830                 for (int j=0 ; j<nbElem; j++)
 831                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s;
 832         return *this;
 833 }
 834
 835 //----------------------------------------------------------------------
 836 const Field&
 837 Field::operator/= ( double s )
 838 //----------------------------------------------------------------------
 839 {
 840         int nbComp=getNumberOfComponents();
 841         int nbElem=getNumberOfElements();
 842         for (int i=0 ; i<nbComp ; i++)
 843                 for (int j=0 ; j<nbElem; j++)
 844                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s;
 845         return *this;
 846 }
 847
 848 //----------------------------------------------------------------------
 849 const Field&
 850 Field::operator-= ( double s )
 851 //----------------------------------------------------------------------
 852 {
 853         int nbComp=getNumberOfComponents();
 854         int nbElem=getNumberOfElements();
 855         for (int i=0 ; i<nbComp ; i++)
 856                 for (int j=0 ; j<nbElem; j++)
 857                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s;
 858         return *this;
 859 }
 860
 861 //----------------------------------------------------------------------
 862 const Field&
 863 Field::operator+= ( double s )
 864 //----------------------------------------------------------------------
 865 {
 866         int nbComp=getNumberOfComponents();
 867         int nbElem=getNumberOfElements();
 868         for (int i=0 ; i<nbComp ; i++)
 869                 for (int j=0 ; j<nbElem; j++)
 870                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s;
 871         return *this;
 872 }
 873
 874 //----------------------------------------------------------------------
 875 void
 876 Field::writeVTK (std::string fileName, bool fromScratch) const
 877 //----------------------------------------------------------------------
 878 {
 879         string fname=fileName+".pvd";
 880         int iter,order;
 881         double time=_field->getTime(iter,order);
 882
 883         if (fromScratch)
 884         {
 885                 ofstream file(fname.c_str()) ;
 886                 file << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\"><Collection>\n" ;
 887                 ostringstream numfile;
 888                 numfile << iter ;
 889                 string filetmp=fileName+"_";
 890                 filetmp=filetmp+numfile.str();
 891                 string ret=_field->writeVTK(filetmp.c_str()) ;
 892                 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
 893                 file << "</Collection></VTKFile>\n" ;
 894                 file.close() ;
 895         }
 896         else
 897         {
 898                 ifstream file1(fname.c_str()) ;
 899                 string contenus;
 900                 getline(file1, contenus, '\0');
 901                 string to_remove="</Collection></VTKFile>";
 902                 size_t m = contenus.find(to_remove);
 903                 size_t n = contenus.find_first_of("\n", m + to_remove.length());
 904                 contenus.erase(m, n - m + 1);
 905                 file1.close() ;
 906                 ofstream file(fname.c_str()) ;
 907                 file << contenus ;
 908                 ostringstream numfile;
 909                 numfile << iter ;
 910                 string filetmp=fileName+"_";
 911                 filetmp=filetmp+numfile.str();
 912                 string ret=_field->writeVTK(filetmp.c_str()) ;
 913                 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
 914                 file << "</Collection></VTKFile>\n" ;
 915                 file.close() ;
 916         }
 917 }
 918
 919 //----------------------------------------------------------------------
 920 void
 921 Field::writeCSV ( const std::string fileName ) const
 922 //----------------------------------------------------------------------
 923 {
 924         int iter,order;
 925         double time=_field->getTime(iter,order);
 926
 927         ostringstream numfile;
 928         numfile << iter ;
 929         string filetmp=fileName+"_";
 930         filetmp=filetmp+numfile.str();
 931         filetmp=filetmp+".csv";
 932         ofstream file(filetmp.c_str()) ;
 933         int dim=_mesh.getSpaceDimension();
 934         int nbElements;
 935         if (getTypeOfField()==CELLS)
 936                 nbElements=_mesh.getNumberOfCells();
 937         else
 938                 nbElements=_mesh.getNumberOfNodes();
 939
 940         if (dim==1)
 941         {
 942                 int nbCompo=getNumberOfComponents();
 943                 if (nbCompo==1)
 944                         file << "x " << _field->getName() << endl;
 945                 else if (nbCompo>1)
 946                 {
 947                         file << "x";
 948                         for (int i=0;i<nbCompo;i++)
 949                                 file << " " << _field->getName() << " Compo " << i+1 << " "<< getInfoOnComponent(i);
 950                         file << endl;
 951                 }
 952                 for (int i=0;i<nbElements;i++)
 953                 {
 954                         if (getTypeOfField()==CELLS)
 955                                 file << _mesh.getCell(i).x() ;
 956                         else
 957                                 file << _mesh.getNode(i).x() ;
 958                         for (int j=0;j<nbCompo;j++)
 959                                 file << " " << getValues()[j+i*nbCompo] ;
 960                         file << endl;
 961                 }
 962         }else if (dim==2)
 963         {
 964                 int nbCompo=getNumberOfComponents();
 965                 if (nbCompo==1)
 966                         file << "x y " << _field->getName() << endl;
 967                 else if (nbCompo>1)
 968                 {
 969                         file << "x y";
 970                         for (int i=0;i<nbCompo;i++)
 971                                 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
 972                         file << endl;
 973                 }
 974                 for (int i=0;i<nbElements;i++)
 975                 {
 976                         if (getTypeOfField()==CELLS)
 977                                 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() ;
 978                         else
 979                                 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() ;
 980                         for (int j=0;j<nbCompo;j++)
 981                                 file << " " << getValues()[j+i*nbCompo] ;
 982                         file << endl;
 983                 }
 984         }else
 985         {
 986                 int nbCompo=getNumberOfComponents();
 987                 if (nbCompo==1)
 988                         file << "x y z " << _field->getName() << endl;
 989                 else if (nbCompo>1)
 990                 {
 991                         file << "x y z";
 992                         for (int i=0;i<nbCompo;i++)
 993                                 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
 994                         file << endl;
 995                 }
 996                 for (int i=0;i<nbElements;i++)
 997                 {
 998                         if (getTypeOfField()==CELLS)
 999                                 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() << " " << _mesh.getCell(i).z();
1000                         else
1001                                 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() << " " << _mesh.getNode(i).z();
1002                         for (int j=0;j<nbCompo;j++)
1003                                 file << " " << getValues()[j+i*nbCompo] ;
1004                         file << endl;
1005                 }
1006         }
1007         file.close() ;
1008 }
1009
1010 //----------------------------------------------------------------------
1011 void
1012 Field::writeMED ( const std::string fileName, bool fromScratch) const
1013 //----------------------------------------------------------------------
1014 {
1015         string fname=fileName+".med";
1016         if (fromScratch)
1017                 MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
1018         else
1019                 MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
1020 }
1021
1022 Field
1023 operator* (double value , const Field& field )
1024 {
1025         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1026         int nbComp=field.getNumberOfComponents();
1027         int nbElem=field.getNumberOfElements();
1028         for (int ielem=0 ; ielem<nbElem; ielem++)
1029                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1030                         fres(ielem, jcomp)=value*field(ielem, jcomp);
1031         return fres;
1032 }
1033
1034 Field
1035 operator* (const Field& field, double value )
1036 {
1037         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1038         int nbComp=field.getNumberOfComponents();
1039         int nbElem=field.getNumberOfElements();
1040         for (int ielem=0 ; ielem<nbElem; ielem++)
1041                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1042                         fres(ielem, jcomp)=value*field(ielem, jcomp);
1043         return fres;
1044 }
1045
1046 Field operator/ (const Field& field, double value)
1047                                 {
1048         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1049         int nbComp=field.getNumberOfComponents();
1050         int nbElem=field.getNumberOfElements();
1051         for (int ielem=0 ; ielem<nbElem; ielem++)
1052                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1053                         fres(ielem, jcomp)=field(ielem, jcomp)/value;
1054         return fres;
1055                                 }
1056
1057 Vector
1058 Field::getValuesOnAllComponents(int elem) const
1059 {
1060         Vector v(getNumberOfComponents());
1061         for(int i=0;i<getNumberOfComponents();i++)
1062                 v[i]=(*this)(elem,i);
1063         return v;
1064 }
1065
1066 Vector
1067 Field::getValuesOnComponent(int compo) const
1068 {
1069         Vector v(getNumberOfElements());
1070         for(int i=0;i<getNumberOfElements();i++)
1071                 v[i]=(*this)(i,compo);
1072         return v;
1073 }
1074
1075 std::ostream& operator<<(std::ostream& out, const Field& field )
1076 {
1077         cout << "Field " << field.getName() << " : " << endl ;
1078         out<< field.getField().retn()->getArray()->repr();
1079         return out;
1080 }