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