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                 DataArrayIdType *desc=DataArrayIdType::New();
  72                 DataArrayIdType *descI=DataArrayIdType::New();
  73                 DataArrayIdType *revDesc=DataArrayIdType::New();
  74                 DataArrayIdType *revDescI=DataArrayIdType::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 Vector
 405 Field::getNormEuclidean() const
 406 {
 407         Vector result(_numberOfComponents);
 408         DoubleTab norm(_numberOfComponents,_field->magnitude()->getArray()->getConstPointer());
 409         result.setValues(norm);
 410
 411         return result;
 412 }
 413
 414 double
 415 Field::max(int component) const
 416 {
 417         if( component >= getNumberOfComponents() || component < 0)
 418                 throw CdmathException("double Field::max() : component number should be between 0 and the field number of components");
 419
 420         double result=-1e100;
 421         for(int i=0; i<getNumberOfElements() ; i++)
 422                 if( result < (*this)(i,component))
 423                         result = (*this)(i,component);
 424
 425         return result;
 426 }
 427
 428 double
 429 Field::min(int component) const
 430 {
 431         if( component >= getNumberOfComponents() || component < 0)
 432                 throw CdmathException("double Field::min() : component number should be between 0 and the field number of components");
 433
 434         double result=1e100;
 435         for(int i=0; i<getNumberOfElements() ; i++)
 436                 if( result > (*this)(i,component))
 437                         result = (*this)(i,component);
 438
 439         return result;
 440 }
 441
 442 string
 443 Field::getInfoOnComponent(int icomp) const
 444 {
 445         return _field->getArray()->getInfoOnComponent(icomp);
 446 }
 447
 448 void
 449 Field::setInfoOnComponent(int icomp, string nameCompo)
 450 {
 451         _field.retn()->getArray()->setInfoOnComponent(icomp,nameCompo);
 452 }
 453
 454 //----------------------------------------------------------------------
 455 Field::Field( const Field & f )
 456 //----------------------------------------------------------------------
 457 {
 458         _mesh=f.getMesh() ;
 459         MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
 460         _field=f1;
 461         _typeField=f.getTypeOfField();
 462 }
 463
 464 //----------------------------------------------------------------------
 465 MCAuto<MEDCouplingFieldDouble>
 466 Field::getField ( void )  const
 467 //----------------------------------------------------------------------
 468 {
 469         return _field ;
 470 }
 471
 472 //----------------------------------------------------------------------
 473 void
 474 Field::setFieldByMEDCouplingFieldDouble ( const MEDCouplingFieldDouble* field )
 475 //----------------------------------------------------------------------
 476 {
 477         MCAuto<MEDCouplingFieldDouble> ff=field->deepCopy();
 478         _field=ff;
 479 }
 480
 481 //----------------------------------------------------------------------
 482 void
 483 Field::setFieldByDataArrayDouble ( const DataArrayDouble* array )
 484 //----------------------------------------------------------------------
 485 {
 486         _field->setArray(const_cast<DataArrayDouble*>(array));
 487 }
 488
 489 //----------------------------------------------------------------------
 490 Vector
 491 Field::integral ( ) const
 492 //----------------------------------------------------------------------
 493 {
 494         int nbComp=_field->getNumberOfComponents();
 495         double res[nbComp];//Pointer containing the inegral of each component
 496         _field->integral(false,res);
 497         Vector result(nbComp);//Vector containing the inegral of each component
 498
 499         for(int i=0; i<nbComp ; i++)
 500                 result(i)=res[i];
 501
 502         return result;
 503 }
 504
 505 //----------------------------------------------------------------------
 506 double
 507 Field::integral ( int compId ) const
 508 //----------------------------------------------------------------------
 509 {
 510         return _field->integral(compId, false);
 511 }
 512
 513 //----------------------------------------------------------------------
 514 Vector
 515 Field::normL1 ( ) const
 516 //----------------------------------------------------------------------
 517 {
 518         int nbComp=_field->getNumberOfComponents();
 519         double res[nbComp];//Pointer containing the L1 norm of each component
 520         _field->normL1(res);
 521         Vector result(nbComp);//Vector containing the L1 norm of each component
 522
 523         for(int i=0; i<nbComp ; i++)
 524                 result(i)=res[i];
 525
 526         return result;
 527 }
 528
 529 //----------------------------------------------------------------------
 530 Vector
 531 Field::normL2 ( ) const
 532 //----------------------------------------------------------------------
 533 {
 534         int nbComp=_field->getNumberOfComponents();
 535         double res[nbComp];//Pointer containing the L2 norm of each component
 536         _field->normL2(res);
 537         Vector result(nbComp);//Vector containing the L2 norm of each component
 538
 539         for(int i=0; i<nbComp ; i++)
 540                 result(i)=res[i];
 541
 542         return result;
 543 }
 544
 545 //----------------------------------------------------------------------
 546 Vector
 547 Field::normMax ( ) const
 548 //----------------------------------------------------------------------
 549 {
 550         int nbComp=_field->getNumberOfComponents();
 551         double res[nbComp];//Pointer containing the L2 norm of each component
 552         _field->normMax(res);
 553         Vector result(nbComp);//Vector containing the L2 norm of each component
 554
 555         for(int i=0; i<nbComp ; i++)
 556                 result(i)=res[i];
 557
 558         return result;
 559 }
 560
 561 Vector
 562 Field::componentMax(Vector & Indices) const
 563 {
 564         int nbComp=_field->getNumberOfComponents();
 565         int nbElems=getNumberOfElements();
 566
 567         Vector result(nbComp);//Vector containing the Linfinity norm of each component
 568
 569         for(int i=0; i<nbElems ; i++)
 570         for(int j=0; j<nbComp ; j++)
 571             if(fabs((*this)(i,j))>result(j))
 572             {
 573                 result(j)=fabs((*this)(i,j));
 574                 Indices(j)=i;
 575             }
 576         return result;
 577 }
 578
 579 //----------------------------------------------------------------------
 580 double&
 581 Field::operator() ( int ielem )
 582 //----------------------------------------------------------------------
 583 {
 584         if(ielem>_field->getNumberOfTuples() || ielem<0)
 585                 throw CdmathException("double& Field::operator(ielem) : ielem>number of values !");
 586         return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
 587 }
 588
 589 //----------------------------------------------------------------------
 590 double&
 591 Field::operator[] ( int ielem )
 592 //----------------------------------------------------------------------
 593 {
 594         if(ielem>_field->getNumberOfTuples() || ielem<0)
 595                 throw CdmathException("double& Field::operator[ielem] : ielem>number of values !");
 596         return _field->getArray()->getPointer()[ielem*_field->getNumberOfComponents()];
 597 }
 598
 599 //----------------------------------------------------------------------
 600 double
 601 Field::operator() ( int ielem ) const
 602 //----------------------------------------------------------------------
 603 {
 604         if(ielem>_field->getNumberOfTuples() || ielem<0)
 605                 throw CdmathException("double Field::operator(ielem) : ielem>number of values !");
 606         return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
 607 }
 608
 609 //----------------------------------------------------------------------
 610 double
 611 Field::operator[] ( int ielem ) const
 612 //----------------------------------------------------------------------
 613 {
 614         if(ielem>_field->getNumberOfTuples() || ielem<0)
 615                 throw CdmathException("double Field::operator[ielem] : ielem>number of values !");
 616         return _field->getArray()->getConstPointer()[ielem*_field->getNumberOfComponents()];
 617 }
 618
 619 //----------------------------------------------------------------------
 620 double&
 621 Field::operator() ( int ielem, int jcomp )
 622 //----------------------------------------------------------------------
 623 {
 624         if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
 625                 throw CdmathException("double& Field::operator( int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
 626         return _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()];
 627 }
 628
 629 //----------------------------------------------------------------------
 630 double
 631 Field::operator() (  int ielem, int jcomp ) const
 632 //----------------------------------------------------------------------
 633 {
 634         if(ielem>_field->getNumberOfTuples() || jcomp>_field->getNumberOfComponents() || ielem<0 || jcomp<0)
 635                 throw CdmathException("double Field::operator(  int ielem, int jcomp ) : ielem>number of values or jcomp>number of components !");
 636         return _field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()];
 637 }
 638
 639 //----------------------------------------------------------------------
 640 void
 641 Field::setTime ( double time, int iter )
 642 //----------------------------------------------------------------------
 643 {
 644         _field->setTime(time,iter,0.0);
 645 }
 646 //----------------------------------------------------------------------
 647 double
 648 Field::getTime ( void ) const
 649 //----------------------------------------------------------------------
 650 {
 651         int a,b;
 652         return _field->getTime(a,b);
 653 }
 654
 655 //----------------------------------------------------------------------
 656 int
 657 Field::getNumberOfElements ( void ) const
 658 //----------------------------------------------------------------------
 659 {
 660         return _field->getNumberOfTuples() ;
 661 }
 662
 663 int
 664 Field::getSpaceDimension( void ) const
 665 {
 666         return _mesh.getSpaceDimension() ;
 667 }
 668
 669 //----------------------------------------------------------------------
 670 int
 671 Field::getNumberOfComponents ( void ) const
 672 //----------------------------------------------------------------------
 673 {
 674         return _field->getNumberOfComponents() ;
 675 }
 676
 677 //----------------------------------------------------------------------
 678 const double*
 679 Field::getValues ( void ) const
 680 //----------------------------------------------------------------------
 681 {
 682         return _field->getArray()->getConstPointer() ;
 683 }
 684
 685 //----------------------------------------------------------------------
 686 const string
 687 Field::getName ( void ) const
 688 //----------------------------------------------------------------------
 689 {
 690         return _field->getName() ;
 691 }
 692
 693 //----------------------------------------------------------------------
 694 const Mesh&
 695 Field::getMesh ( void ) const
 696 //----------------------------------------------------------------------
 697 {
 698         return _mesh ;
 699 }
 700
 701 //----------------------------------------------------------------------
 702 EntityType
 703 Field::getTypeOfField ( void ) const
 704 //----------------------------------------------------------------------
 705 {
 706         return _typeField;
 707 }
 708
 709 //----------------------------------------------------------------------
 710 void
 711 Field::setName ( const string fieldName )
 712 //----------------------------------------------------------------------
 713 {
 714         _field->setName(fieldName.c_str()) ;
 715 }
 716
 717
 718 //----------------------------------------------------------------------
 719 Field
 720 Field::operator+ ( const Field& f ) const
 721 //----------------------------------------------------------------------
 722 {
 723         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 724         //throw CdmathException("Field::operator+ : Field addition requires identical meshes");
 725         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 726         if(f.getTypeOfField() != getTypeOfField())
 727                 throw CdmathException("Field::operator+ : Field addition requires identical field types (CELLS, NODES or FACES");
 728         if(f.getNumberOfComponents() != getNumberOfComponents())
 729                 throw CdmathException("Field::operator+ : Field addition requires identical number of components");
 730
 731         Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
 732         int nbComp=f.getNumberOfComponents();
 733         int nbElem=f.getNumberOfElements();
 734         for (int ielem=0 ; ielem<nbElem; ielem++)
 735                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
 736                         fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp);
 737         return fres;
 738 }
 739
 740 //----------------------------------------------------------------------
 741 Field
 742 Field::operator- ( const Field& f ) const
 743 //----------------------------------------------------------------------
 744 {
 745         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 746         //throw CdmathException("Field::operator- : Field subtraction requires identical meshes");
 747         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 748         if(f.getTypeOfField() != getTypeOfField())
 749                 throw CdmathException("Field::operator- : Field subtraction requires identical field types (CELLS, NODES or FACES");
 750         if(f.getNumberOfComponents() != getNumberOfComponents())
 751                 throw CdmathException("Field::operator- : Field subtraction requires identical number of components");
 752
 753         Field fres(getName(),f.getTypeOfField(),f.getMesh(),f.getNumberOfComponents(),f.getTime());
 754         int nbComp=f.getNumberOfComponents();
 755         int nbElem=f.getNumberOfElements();
 756         for (int ielem=0 ; ielem<nbElem; ielem++)
 757                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
 758                         fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp);
 759         return fres;
 760 }
 761
 762 //----------------------------------------------------------------------
 763 const Field&
 764 Field::operator= ( const Field& f )
 765 //----------------------------------------------------------------------
 766 {
 767         _mesh=f.getMesh() ;
 768         _typeField=f.getTypeOfField() ;
 769         _numberOfComponents=f.getNumberOfComponents();
 770         _time=f.getTime();
 771         _fieldName=f.getName();
 772         MCAuto<MEDCouplingFieldDouble> f1=f.getField()->deepCopy();
 773         _field=f1;
 774         return *this;
 775 }
 776
 777 Field Field::deepCopy( ) const
 778 {
 779     Field F(getName(), getTypeOfField(), getMesh(), getNumberOfComponents(), getTime()) ;
 780         MCAuto<MEDCouplingFieldDouble> f1=getField()->deepCopy();
 781         F.setFieldByMEDCouplingFieldDouble(f1);
 782
 783     return F;
 784 }
 785
 786
 787 //----------------------------------------------------------------------
 788 const Field&
 789 Field::operator+= ( const Field& f )
 790 //----------------------------------------------------------------------
 791 {
 792         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 793         //throw CdmathException("Field::operator+= : Field addition requires identical meshes");
 794         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 795         if(f.getTypeOfField() != getTypeOfField())
 796                 throw CdmathException("Field::operator+= : Field addition requires identical field types (CELLS, NODES or FACES");
 797         if(f.getNumberOfComponents() != getNumberOfComponents())
 798                 throw CdmathException("Field::operator+= : Field addition requires identical number of components");
 799
 800         _field->setMesh(f.getField()->getMesh());
 801         (*_field)+=(*f.getField());
 802         return *this;
 803 }
 804
 805 //----------------------------------------------------------------------
 806 const Field&
 807 Field::operator-= ( const Field& f )
 808 //----------------------------------------------------------------------
 809 {
 810         //if(f.getMesh().getMEDCouplingMesh() != _mesh.getMEDCouplingMesh())
 811         //throw CdmathException("Field::operator-= : Field subtraction requires identical meshes");
 812         _mesh.getMEDCouplingMesh()->checkFastEquivalWith(f.getMesh().getMEDCouplingMesh(),1e-6);
 813         if(f.getTypeOfField() != getTypeOfField())
 814                 throw CdmathException("Field::operator-= : Field subtraction requires identical field types (CELLS, NODES or FACES");
 815         if(f.getNumberOfComponents() != getNumberOfComponents())
 816                 throw CdmathException("Field::operator-= : Field subtraction requires identical number of components");
 817
 818         _field->setMesh(f.getField()->getMesh());
 819         (*_field)-=(*f.getField());
 820         return *this;
 821 }
 822
 823 //----------------------------------------------------------------------
 824 const Field&
 825 Field::operator*= ( double s )
 826 //----------------------------------------------------------------------
 827 {
 828         int nbComp=getNumberOfComponents();
 829         int nbElem=getNumberOfElements();
 830         for (int i=0 ; i<nbComp ; i++)
 831                 for (int j=0 ; j<nbElem; j++)
 832                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s;
 833         return *this;
 834 }
 835
 836 //----------------------------------------------------------------------
 837 const Field&
 838 Field::operator/= ( double s )
 839 //----------------------------------------------------------------------
 840 {
 841         int nbComp=getNumberOfComponents();
 842         int nbElem=getNumberOfElements();
 843         for (int i=0 ; i<nbComp ; i++)
 844                 for (int j=0 ; j<nbElem; j++)
 845                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s;
 846         return *this;
 847 }
 848
 849 //----------------------------------------------------------------------
 850 const Field&
 851 Field::operator-= ( double s )
 852 //----------------------------------------------------------------------
 853 {
 854         int nbComp=getNumberOfComponents();
 855         int nbElem=getNumberOfElements();
 856         for (int i=0 ; i<nbComp ; i++)
 857                 for (int j=0 ; j<nbElem; j++)
 858                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s;
 859         return *this;
 860 }
 861
 862 //----------------------------------------------------------------------
 863 const Field&
 864 Field::operator+= ( double s )
 865 //----------------------------------------------------------------------
 866 {
 867         int nbComp=getNumberOfComponents();
 868         int nbElem=getNumberOfElements();
 869         for (int i=0 ; i<nbComp ; i++)
 870                 for (int j=0 ; j<nbElem; j++)
 871                         _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s;
 872         return *this;
 873 }
 874
 875 //----------------------------------------------------------------------
 876 void
 877 Field::writeVTK (std::string fileName, bool fromScratch) const
 878 //----------------------------------------------------------------------
 879 {
 880         string fname=fileName+".pvd";
 881         int iter,order;
 882         double time=_field->getTime(iter,order);
 883
 884         if (fromScratch)
 885         {
 886                 ofstream file(fname.c_str()) ;
 887                 file << "<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\"><Collection>\n" ;
 888                 ostringstream numfile;
 889                 numfile << iter ;
 890                 string filetmp=fileName+"_";
 891                 filetmp=filetmp+numfile.str();
 892                 string ret=_field->writeVTK(filetmp.c_str()) ;
 893                 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
 894                 file << "</Collection></VTKFile>\n" ;
 895                 file.close() ;
 896         }
 897         else
 898         {
 899                 ifstream file1(fname.c_str()) ;
 900                 string contenus;
 901                 getline(file1, contenus, '\0');
 902                 string to_remove="</Collection></VTKFile>";
 903                 size_t m = contenus.find(to_remove);
 904                 size_t n = contenus.find_first_of("\n", m + to_remove.length());
 905                 contenus.erase(m, n - m + 1);
 906                 file1.close() ;
 907                 ofstream file(fname.c_str()) ;
 908                 file << contenus ;
 909                 ostringstream numfile;
 910                 numfile << iter ;
 911                 string filetmp=fileName+"_";
 912                 filetmp=filetmp+numfile.str();
 913                 string ret=_field->writeVTK(filetmp.c_str()) ;
 914                 file << "<DataSet timestep=\""<< time << "\" group=\"\" part=\"0\" file=\"" << ret << "\"/>\n" ;
 915                 file << "</Collection></VTKFile>\n" ;
 916                 file.close() ;
 917         }
 918 }
 919
 920 //----------------------------------------------------------------------
 921 void
 922 Field::writeCSV ( const std::string fileName ) const
 923 //----------------------------------------------------------------------
 924 {
 925         int iter,order;
 926         double time=_field->getTime(iter,order);
 927
 928         ostringstream numfile;
 929         numfile << iter ;
 930         string filetmp=fileName+"_";
 931         filetmp=filetmp+numfile.str();
 932         filetmp=filetmp+".csv";
 933         ofstream file(filetmp.c_str()) ;
 934         int dim=_mesh.getSpaceDimension();
 935         int nbElements;
 936         if (getTypeOfField()==CELLS)
 937                 nbElements=_mesh.getNumberOfCells();
 938         else
 939                 nbElements=_mesh.getNumberOfNodes();
 940
 941         if (dim==1)
 942         {
 943                 int nbCompo=getNumberOfComponents();
 944                 if (nbCompo==1)
 945                         file << "x " << _field->getName() << endl;
 946                 else if (nbCompo>1)
 947                 {
 948                         file << "x";
 949                         for (int i=0;i<nbCompo;i++)
 950                                 file << " " << _field->getName() << " Compo " << i+1 << " "<< getInfoOnComponent(i);
 951                         file << endl;
 952                 }
 953                 for (int i=0;i<nbElements;i++)
 954                 {
 955                         if (getTypeOfField()==CELLS)
 956                                 file << _mesh.getCell(i).x() ;
 957                         else
 958                                 file << _mesh.getNode(i).x() ;
 959                         for (int j=0;j<nbCompo;j++)
 960                                 file << " " << getValues()[j+i*nbCompo] ;
 961                         file << endl;
 962                 }
 963         }else if (dim==2)
 964         {
 965                 int nbCompo=getNumberOfComponents();
 966                 if (nbCompo==1)
 967                         file << "x y " << _field->getName() << endl;
 968                 else if (nbCompo>1)
 969                 {
 970                         file << "x y";
 971                         for (int i=0;i<nbCompo;i++)
 972                                 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
 973                         file << endl;
 974                 }
 975                 for (int i=0;i<nbElements;i++)
 976                 {
 977                         if (getTypeOfField()==CELLS)
 978                                 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() ;
 979                         else
 980                                 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() ;
 981                         for (int j=0;j<nbCompo;j++)
 982                                 file << " " << getValues()[j+i*nbCompo] ;
 983                         file << endl;
 984                 }
 985         }else
 986         {
 987                 int nbCompo=getNumberOfComponents();
 988                 if (nbCompo==1)
 989                         file << "x y z " << _field->getName() << endl;
 990                 else if (nbCompo>1)
 991                 {
 992                         file << "x y z";
 993                         for (int i=0;i<nbCompo;i++)
 994                                 file << " " << _field->getName() << " Compo " << i+1<< " "<< getInfoOnComponent(i);
 995                         file << endl;
 996                 }
 997                 for (int i=0;i<nbElements;i++)
 998                 {
 999                         if (getTypeOfField()==CELLS)
1000                                 file << _mesh.getCell(i).x() << " " << _mesh.getCell(i).y() << " " << _mesh.getCell(i).z();
1001                         else
1002                                 file << _mesh.getNode(i).x() << " " << _mesh.getNode(i).y() << " " << _mesh.getNode(i).z();
1003                         for (int j=0;j<nbCompo;j++)
1004                                 file << " " << getValues()[j+i*nbCompo] ;
1005                         file << endl;
1006                 }
1007         }
1008         file.close() ;
1009 }
1010
1011 //----------------------------------------------------------------------
1012 void
1013 Field::writeMED ( const std::string fileName, bool fromScratch) const
1014 //----------------------------------------------------------------------
1015 {
1016         string fname=fileName+".med";
1017         if (fromScratch)
1018                 MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
1019         else
1020                 MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
1021 }
1022
1023 Field
1024 operator* (double value , const Field& field )
1025 {
1026         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1027         int nbComp=field.getNumberOfComponents();
1028         int nbElem=field.getNumberOfElements();
1029         for (int ielem=0 ; ielem<nbElem; ielem++)
1030                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1031                         fres(ielem, jcomp)=value*field(ielem, jcomp);
1032         return fres;
1033 }
1034
1035 Field
1036 operator* (const Field& field, double value )
1037 {
1038         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1039         int nbComp=field.getNumberOfComponents();
1040         int nbElem=field.getNumberOfElements();
1041         for (int ielem=0 ; ielem<nbElem; ielem++)
1042                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1043                         fres(ielem, jcomp)=value*field(ielem, jcomp);
1044         return fres;
1045 }
1046
1047 Field operator/ (const Field& field, double value)
1048                                 {
1049         Field fres(field.getName(),field.getTypeOfField(),field.getMesh(),field.getNumberOfComponents(),field.getTime());
1050         int nbComp=field.getNumberOfComponents();
1051         int nbElem=field.getNumberOfElements();
1052         for (int ielem=0 ; ielem<nbElem; ielem++)
1053                 for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
1054                         fres(ielem, jcomp)=field(ielem, jcomp)/value;
1055         return fres;
1056                                 }
1057
1058 Vector
1059 Field::getValuesOnAllComponents(int elem) const
1060 {
1061         Vector v(getNumberOfComponents());
1062         for(int i=0;i<getNumberOfComponents();i++)
1063                 v[i]=(*this)(elem,i);
1064         return v;
1065 }
1066
1067 Vector
1068 Field::getValuesOnComponent(int compo) const
1069 {
1070         Vector v(getNumberOfElements());
1071         for(int i=0;i<getNumberOfElements();i++)
1072                 v[i]=(*this)(i,compo);
1073         return v;
1074 }
1075
1076 std::ostream& operator<<(std::ostream& out, const Field& field )
1077 {
1078         cout << "Field " << field.getName() << " : " << endl ;
1079         out<< field.getField().retn()->getArray()->repr();
1080         return out;
1081 }