Start using profile to deal with 1) boundary conditions and 2) mesh deletion in memory

diff --git a/CDMATH/mesh/inc/Field.hxx b/CDMATH/mesh/inc/Field.hxx
index 542a68828597d9764bd0e6f49bf72ce216199ce1..489588c16c57b5cb8d2f0d96eab25d694082e7c6 100755 (executable)
--- a/CDMATH/mesh/inc/Field.hxx
+++ b/CDMATH/mesh/inc/Field.hxx
@@ -12,6 +12,7 @@ namespace MEDCoupling
 {
   class MEDCouplingFieldDouble;
   class DataArrayDouble;
+  class MEDFileField1TS;
 }
 
 #include "DoubleTab.hxx"
@@ -72,8 +73,7 @@ class Field
      */
     Field( const std::string filename, EntityType fieldType,
            const std::string & fieldName = "",
-           int iteration = -1, int order = -1, int meshLevel=0,
-           int numberOfComponents=1, double time=0.0);
+           int iteration = -1, int order = -1, int meshLevel=0);
   
     /**
      * constructor with data
@@ -251,6 +251,16 @@ class Field
 
     void setFieldByDataArrayDouble ( const MEDCoupling::DataArrayDouble* array );
 
+    /** 
+        * \brief Delete the medcoupling mesh to save memory space
+     */
+    void deleteMEDCouplingUMesh();
+    
+    /** 
+     * \brief Returns true iff an unstructured mesh has been loaded
+     */
+     bool meshNotDeleted() const {return _mesh.meshNotDeleted();}
+    
     Vector getNormEuclidean( void ) const ;
 
     double max( int component=0 ) const ;
@@ -357,6 +367,7 @@ class Field
        int _numberOfComponents;
        double _time;
        std::string _fieldName;
+       MEDCoupling::MEDFileField1TS *_ff ;//To save the field when the mesh has been deleted
        
     private:
 

diff --git a/CDMATH/mesh/src/Field.cxx b/CDMATH/mesh/src/Field.cxx
index 268af9766c7547116c41330d5ada081e9e6bd836..3252115340379bb654ac2d38b01a44df313cdd5f 100755 (executable)
--- a/CDMATH/mesh/src/Field.cxx
+++ b/CDMATH/mesh/src/Field.cxx
@@ -10,6 +10,7 @@
 #include "Face.hxx"
 #include "Field.hxx"
 #include "MEDFileMesh.hxx"
+#include "MEDFileField1TS.hxx"
 #include "MEDLoader.hxx"
 #include "MEDCouplingUMesh.hxx"
 #include "MEDCouplingFieldDouble.hxx"
@@ -29,6 +30,7 @@ Field::Field( EntityType typeField )
 //----------------------------------------------------------------------
 {
        _field=NULL;
+       _ff=NULL;
        _typeField=typeField;
        _numberOfComponents=0;
 }
@@ -45,6 +47,7 @@ Field::~Field( void )
 Field::Field(const std::string fieldName, EntityType type, const Mesh& mesh, int numberOfComponents, double time)
 {
        _field = NULL;
+       _ff=NULL;
        _mesh=Mesh(mesh);
        _typeField=type;
        _numberOfComponents=numberOfComponents;
@@ -89,19 +92,16 @@ void Field::buildFieldMemoryStructure()
        _field->setArray(array);
        _field->setTime(_time,0,0);
        array->decrRef();
-       mu->decrRef();
 }
 
 Field::Field( const std::string filename, EntityType type,
                const std::string & fieldName,
-               int iteration, int order, int meshLevel,
-               int numberOfComponents, double time)
+               int iteration, int order, int meshLevel)
 {
        _field = NULL;
+       _ff=NULL;
        _mesh=Mesh(filename + ".med", meshLevel);
        _typeField=type;
-       _numberOfComponents=numberOfComponents;
-       _time=time;
        _fieldName=fieldName;
 
        readFieldMed(filename, type, fieldName, iteration, order);
@@ -111,6 +111,7 @@ Field::Field(const std::string meshFileName, EntityType type, const std::vector<
                const std::string & fieldName, int meshLevel, double time )
 {
        _field = NULL;
+       _ff=NULL;
        _mesh=Mesh(meshFileName + ".med", meshLevel);
        _typeField=type;
        _numberOfComponents=Vconstant.size();
@@ -124,11 +125,12 @@ Field::Field(const std::string meshFileName, EntityType type, const std::vector<
 
        for (int ielem=0 ; ielem<nbelem; ielem++)
                for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
-                       _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
+                       _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp];
 }
 Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::string & fieldName, double time)
 {
        _field = NULL;
+       _ff=NULL;
        _mesh=Mesh(M);
        _typeField=type;
        _numberOfComponents=Vconstant.size();
@@ -142,11 +144,12 @@ Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::
 
        for (int ielem=0 ; ielem<nbelem; ielem++)
                for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
-                       _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
+                       _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp];
 }
 Field::Field(const Mesh& M, EntityType type, const vector<double> Vconstant, const std::string & fieldName, double time) 
 {
        _field = NULL;
+       _ff=NULL;
        _mesh=Mesh(M);
        _typeField=type;
        _numberOfComponents=Vconstant.size();
@@ -160,7 +163,7 @@ Field::Field(const Mesh& M, EntityType type, const vector<double> Vconstant, con
 
        for (int ielem=0 ; ielem<nbelem; ielem++)
                for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
-                       _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
+                       _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp];
 }
 Field::Field( int nDim, const vector<double> Vconstant, EntityType type, 
                double xmin, double xmax, int nx, string leftSide, string rightSide,
@@ -169,6 +172,7 @@ Field::Field( int nDim, const vector<double> Vconstant, EntityType type,
                const std::string & fieldName, double time,double epsilon)
 {
        _field = NULL;
+       _ff=NULL;
        _typeField=type;
        _numberOfComponents=Vconstant.size();
        _time=time;
@@ -206,7 +210,7 @@ Field::Field( int nDim, const vector<double> Vconstant, EntityType type,
 
        for (int ielem=0 ; ielem<nbelem; ielem++)
                for (int jcomp=0 ; jcomp<nbcomp ; jcomp++)
-                       _field->getArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp];
+                       _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp];
 }
 Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos,
                EntityType type, int direction, const std::string & fieldName, double time)
@@ -215,6 +219,7 @@ Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double d
                throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
 
        _field = NULL;
+       _ff=NULL;
        _mesh=Mesh(M);
        _typeField=type;
        _numberOfComponents=VV_Left.getNumberOfRows();
@@ -240,9 +245,9 @@ Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double d
 
                for (int i=0; i< nbcomp; i++)
                        if (component_value< disc_pos )
-                               _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Left[i];
+                               _field->getArray()->getPointer()[j+i*nbcomp] = VV_Left[i];
                        else
-                               _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Right[i];
+                               _field->getArray()->getPointer()[j+i*nbcomp] = VV_Right[i];
        }
 }
 Field::Field( int nDim, const vector<double> VV_Left, vector<double> VV_Right, 
@@ -292,6 +297,7 @@ Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius,
        }
 
        _field = NULL;
+       _ff=NULL;
        _mesh=Mesh(M);
        _typeField=type;
        _numberOfComponents=Vout.size();
@@ -318,10 +324,10 @@ Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius,
                }
                if((currentPoint-Center).norm()<radius)
                        for(int j=0;j<nbcomp;j++)
-                               _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vin[j];
+                               _field->getArray()->getPointer()[j+i*nbcomp]=Vin[j];
                else
                        for(int j=0;j<nbcomp;j++)
-                               _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vout[j];
+                               _field->getArray()->getPointer()[j+i*nbcomp]=Vout[j];
        }
 }
 
@@ -344,23 +350,33 @@ Field::readFieldMed( const std::string & fileNameRadical,
        size_t iField = 0;
        std::string attributedFieldName;
        _field = NULL;
+       _ff=NULL;
 
        // Get the name of the right field that we will attribute to the Field.
        if (fieldName == "") {
                if (fieldNames.size() > 0)
+               {
+                       cout<<"Warning : No field name imposed, taking the first field name found : "<< fieldNames[0]<<" in file "<<completeFileName<<endl;;
                        attributedFieldName = fieldNames[0];
-               else {
+               }
+               else 
+               {
                        std::ostringstream message;
-                       message << "No field in file " << completeFileName;
+                       message << "Warning : No field found in file " << completeFileName;
                        throw CdmathException(message.str().c_str());
                }
        }
-       else {
-               for (; iField < fieldNames.size(); iField++)
-                       if (fieldName == fieldNames[iField]) break;
-
-               if (iField < fieldNames.size())
+       else 
+       {
+               int mycount = std::count(fieldNames.begin(), fieldNames.end(), fieldName);
+               
+               if( mycount>0 )
+               {
                        attributedFieldName = fieldName;
+
+                       if( mycount> 1 )
+                               cout<<"Warning : " << mycount << " fields are associated to the name " << fieldName <<" in file "<<completeFileName<<endl;
+               }
                else {
                        std::ostringstream message;
                        message << "No field named " << fieldName << " in file " << completeFileName;
@@ -369,38 +385,52 @@ Field::readFieldMed( const std::string & fileNameRadical,
        }
 
        // Get the name of the right mesh that we will attribute to the Field.
-       std::vector<std::string> meshNames
-       = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName);
+       std::vector<std::string> meshNames      = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName);
        if (meshNames.size() == 0) {
                std::ostringstream message;
-               message << "No mesh associated to " << fieldName
-                               << " in file " << completeFileName;
+               message << "No mesh associated to " << fieldName<< " in file " << completeFileName;
                throw CdmathException(message.str().c_str());
        }
+       else
+               if( meshNames.size() > 1 )
+               {
+                       cout<<"Warning : " << meshNames.size() << " meshes are associated to field named " << fieldName <<" in file "<<completeFileName<<endl;
+                       cout<<"Mesh names are : ";
+                       for (int iMesh=0; iMesh < meshNames.size(); iMesh++)
+                               cout<< meshNames[iMesh]<<", ";
+                       cout<<endl<<"Taking the mesh with name "<< meshNames[0]<<endl;                  
+               }
+
        std::string attributedMeshName = meshNames[0];
 
        // Create Field.
-       MEDCoupling::TypeOfField medFieldType[3] = { ON_CELLS, ON_NODES, ON_CELLS };
-       switch (type) {
-       case CELLS:
-               _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( 
-                                       MEDCoupling::ReadFieldCell( completeFileName,
-                                       attributedMeshName, 0,
-                                       attributedFieldName, iteration, order) );
-               break;
-       case NODES:
-               _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
-                                       MEDCoupling::ReadFieldNode( completeFileName,
-                                       attributedMeshName, 0,
-                                       attributedFieldName, iteration, order) );
-               break;
-       case FACES:
-               _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( 
-                                       MEDCoupling::ReadFieldCell( completeFileName,
-                                       attributedMeshName, -1,
-                                       attributedFieldName, iteration, order) );
-               break;
+       switch (type) 
+       {
+               case CELLS:
+                       _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( 
+                                               MEDCoupling::ReadFieldCell( completeFileName,
+                                               attributedMeshName, 0,
+                                               attributedFieldName, iteration, order) );
+                       break;
+               case NODES:
+                       _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > (
+                                               MEDCoupling::ReadFieldNode( completeFileName,
+                                               attributedMeshName, 0,
+                                               attributedFieldName, iteration, order) );
+                       break;
+               case FACES:
+                       _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( 
+                                               MEDCoupling::ReadFieldCell( completeFileName,
+                                               attributedMeshName, -1,
+                                               attributedFieldName, iteration, order) );
+                       break;
        }
+
+       //Read and store the number of components
+       _numberOfComponents = _field->getNumberOfComponents() ;
+       _time = _field->getTime(iteration, order);
+
+       cout<<"Found field " << fieldName << " in file " << completeFileName <<endl;
 }
 
 
@@ -818,7 +848,7 @@ Field::operator+ ( const Field& f ) const
        int nbElem=f.getNumberOfElements();
        for (int ielem=0 ; ielem<nbElem; ielem++)
                for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
-                       fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp);
+                       fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*nbComp]+f(ielem, jcomp);
        return fres;
 }
 
@@ -840,7 +870,7 @@ Field::operator- ( const Field& f ) const
        int nbElem=f.getNumberOfElements();
        for (int ielem=0 ; ielem<nbElem; ielem++)
                for (int jcomp=0 ; jcomp<nbComp ; jcomp++)
-                       fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp);
+                       fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*nbComp]-f(ielem, jcomp);
        return fres;
 }
 
@@ -914,7 +944,7 @@ Field::operator*= ( double s )
        int nbElem=getNumberOfElements();
        for (int i=0 ; i<nbComp ; i++)
                for (int j=0 ; j<nbElem; j++)
-                       _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s;
+                       _field->getArray()->getPointer()[i+j*nbComp]*=s;
        return *this;
 }
 
@@ -927,7 +957,7 @@ Field::operator/= ( double s )
        int nbElem=getNumberOfElements();
        for (int i=0 ; i<nbComp ; i++)
                for (int j=0 ; j<nbElem; j++)
-                       _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s;
+                       _field->getArray()->getPointer()[i+j*nbComp]/=s;
        return *this;
 }
 
@@ -940,7 +970,7 @@ Field::operator-= ( double s )
        int nbElem=getNumberOfElements();
        for (int i=0 ; i<nbComp ; i++)
                for (int j=0 ; j<nbElem; j++)
-                       _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s;
+                       _field->getArray()->getPointer()[i+j*nbComp]-=s;
        return *this;
 }
 
@@ -953,7 +983,7 @@ Field::operator+= ( double s )
        int nbElem=getNumberOfElements();
        for (int i=0 ; i<nbComp ; i++)
                for (int j=0 ; j<nbElem; j++)
-                       _field->getArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s;
+                       _field->getArray()->getPointer()[i+j*nbComp]+=s;
        return *this;
 }
 
@@ -962,6 +992,9 @@ void
 Field::writeVTK (std::string fileName, bool fromScratch) const
 //----------------------------------------------------------------------
 {
+       if( !_mesh.isStructured() && !_mesh.meshNotDeleted() )
+               throw CdmathException("Field::writeVTK : Cannot save field in VTK format : unstructured mesh with no MEDCouplingUMesh loaded. Use med format.");
+
        string fname=fileName+".pvd";
        int iter,order;
        double time=_field->getTime(iter,order);
@@ -1099,10 +1132,27 @@ Field::writeMED ( const std::string fileName, bool fromScratch) const
 //----------------------------------------------------------------------
 {
        string fname=fileName+".med";
-       if (fromScratch)
-               MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
-       else
-               MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
+       
+       if(_mesh.isStructured() || _mesh.meshNotDeleted())
+               if (fromScratch)
+                       MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
+               else
+                       MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
+       else//The mesh has ben deleted, use _ff instead of _field to save the values
+       {
+               //MEDFileUMesh * meshMEDFile = MEDFileUMesh::New();
+               //meshMEDFile->setMeshAtLevel(0,_field->getMesh()->buildUnstructured());
+               //meshMEDFile->write(fname.c_str(), fromScratch);
+               //MEDCoupling::WriteUMesh(fname.c_str(),_field->getMesh()->buildUnstructured(),fromScratch);
+               //MEDCoupling::WriteMesh(fname.c_str(),_field->getMesh(),fromScratch);
+           //MEDCoupling::MEDCouplingUMesh* fmesh = dynamic_cast<MEDCoupling::MEDCouplingUMesh*> (_field->getMesh()->deepCopy());
+               //cout<<" checkConsecutiveCellTypes : "<< fmesh->checkConsecutiveCellTypes() <<endl;
+               //cout<<" advancedRepr() : "<< fmesh->advancedRepr() <<endl;
+               //cout<<" checkConsecutiveCellTypes : "<< _field->getMesh()->buildUnstructured()->checkConsecutiveCellTypes()<<endl;
+               MEDFileField1TS *ff=MEDFileField1TS::New();
+               _ff->setFieldNoProfileSBT(  _field );
+               _ff->write(fname.c_str(), fromScratch);
+       }
 }
 
 Field
@@ -1173,3 +1223,8 @@ std::ostream& operator<<(std::ostream& out, const Field& field )
        out<< field.getField().retn()->getArray()->repr();
        return out;
 }
+
+void Field::deleteMEDCouplingUMesh()
+{ 
+       return _mesh.deleteMEDCouplingUMesh();
+}