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