Better handling of field support (CELLS vs NODES)

diff --git a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx
index 7810a39d3cd48eaf8d428abedeeb8d8d6ae99d48..7bd065ea5841f484e7acc18fc3ebbc9db3aa66a1 100755 (executable)
--- a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx
+++ b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx
@@ -217,7 +217,7 @@ protected :
     std::vector< int > _dirichletNodeIds;/* List of boundary nodes with Dirichlet BC */
 
     /*********** Functions for finite element method ***********/
-    Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
+    static Vector gradientNodal(Matrix M, vector< double > v);//gradient of nodal shape functions
        double computeDiffusionMatrixFE(bool & stop);
     static int fact(int n);
     static int unknownNodeIndex(int globalIndex, std::vector< int > dirichletNodes);

diff --git a/CoreFlows/Models/src/ProblemCoreFlows.cxx b/CoreFlows/Models/src/ProblemCoreFlows.cxx
index c19cc8d9c4613e300870a5168573c11ac630127e..aea9a503b663c73df70cec150b87bf718b0a881e 100755 (executable)
--- a/CoreFlows/Models/src/ProblemCoreFlows.cxx
+++ b/CoreFlows/Models/src/ProblemCoreFlows.cxx
@@ -124,6 +124,12 @@ void ProblemCoreFlows::setInitialField(const Field &VV)
                _runLogFile->close();
                throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect number of components");
        }
+       if(_FECalculation && VV.getTypeOfField()!=NODES)
+       {
+               *_runLogFile<<"ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Element method"<<endl;
+               _runLogFile->close();
+               throw CdmathException("ProblemCoreFlows::setInitialField: Initial field has incorrect support : should be on nodes for the Finite Element method");
+       }
 
        _VV=VV;
        _VV.setName("SOLVERLAB results");
@@ -622,10 +628,7 @@ ProblemCoreFlows::getEigenvectorsField(int nev, EPSWhich which, double tol) cons
   MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
   Field my_eigenfield;
   
-  if(_FECalculation)
-    my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
-  else
-    my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
+  my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
 
   my_eigenfield.setFieldByDataArrayDouble(d);
   

diff --git a/CoreFlows/Models/src/ProblemFluid.cxx b/CoreFlows/Models/src/ProblemFluid.cxx
index 8f497d37752851cda25c3a8a8b599e1b2fa626c5..731b5a72c4c30b159f06af7d3d87a86e9f21b0a4 100755 (executable)
--- a/CoreFlows/Models/src/ProblemFluid.cxx
+++ b/CoreFlows/Models/src/ProblemFluid.cxx
@@ -42,11 +42,18 @@ void ProblemFluid::setNumericalScheme(SpaceScheme spaceScheme, TimeScheme timeSc
 
 void ProblemFluid::initialize()
 {
-       if(!_initialDataSet){
+       if(!_initialDataSet)
+       {
                *_runLogFile<<"ProblemFluid::initialize() set initial data first"<<endl;
                _runLogFile->close();
                throw CdmathException("ProblemFluid::initialize() set initial data first");
        }
+       else if (_VV.getTypeOfField() != CELLS)
+       {
+               *_runLogFile<<"Initial data should be a field on CELLS, not NODES, neither FACES"<<endl;
+               _runLogFile->close();
+               throw CdmathException("ProblemFluid::initialize() Initial data should be a field on CELLS, not NODES, neither FACES");
+       }
        cout << "Number of Phases = " << _nbPhases << " mesh dimension = "<<_Ndim<<" number of variables = "<<_nVar<<endl;
        *_runLogFile << "Number of Phases = " << _nbPhases << " spaceDim= "<<_Ndim<<" number of variables= "<<_nVar<<endl;
 

diff --git a/CoreFlows/Models/src/StationaryDiffusionEquation.cxx b/CoreFlows/Models/src/StationaryDiffusionEquation.cxx
index a4661a3115da9d6047f1e0979283d6331757fe63..44ae333bfbfbb811805c547e9eb6f3c044f18d95 100755 (executable)
--- a/CoreFlows/Models/src/StationaryDiffusionEquation.cxx
+++ b/CoreFlows/Models/src/StationaryDiffusionEquation.cxx
@@ -151,29 +151,15 @@ void StationaryDiffusionEquation::initialize()
        /**************** Field creation *********************/
 
        if(!_heatPowerFieldSet){
-        if(_FECalculation){
-            _heatPowerField=Field("Heat power",NODES,_mesh,1);
-            for(int i =0; i<_Nnodes; i++)
-                _heatPowerField(i) = _heatSource;
-        }
-        else{
-            _heatPowerField=Field("Heat power",CELLS,_mesh,1);
-            for(int i =0; i<_Nmailles; i++)
-                _heatPowerField(i) = _heatSource;
-        }
+               _heatPowerField=Field("Heat power",_VV.getTypeOfField(),_mesh,1);
+               for(int i =0; i<_VV.getNumberOfElements(); i++)
+                       _heatPowerField(i) = _heatSource;
         _heatPowerFieldSet=true;
     }
        if(!_fluidTemperatureFieldSet){
-        if(_FECalculation){
-            _fluidTemperatureField=Field("Fluid temperature",NODES,_mesh,1);
-            for(int i =0; i<_Nnodes; i++)
-                _fluidTemperatureField(i) = _fluidTemperature;
-        }
-        else{
-            _fluidTemperatureField=Field("Fluid temperature",CELLS,_mesh,1);
-            for(int i =0; i<_Nmailles; i++)
-                _fluidTemperatureField(i) = _fluidTemperature;
-        }
+               _fluidTemperatureField=Field("Fluid temperature",_VV.getTypeOfField(),_mesh,1);
+               for(int i =0; i<_VV.getNumberOfElements(); i++)
+                       _fluidTemperatureField(i) = _fluidTemperature;
         _fluidTemperatureFieldSet=true;
        }
 
@@ -293,17 +279,22 @@ double StationaryDiffusionEquation::computeTimeStep(bool & stop){
        return _dt_src;
 }
 
-Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values){
-    vector< Matrix > matrices(_Ndim);
+Vector StationaryDiffusionEquation::gradientNodal(Matrix M, vector< double > values)
+{
+       if(! M.isSquare() )
+               throw CdmathException("DiffusionEquation::gradientNodal Matrix M should be square !!!");
+               
+       int Ndim = M.getNumberOfRows();
+    vector< Matrix > matrices(Ndim);
     
-    for (int idim=0; idim<_Ndim;idim++){
+    for (int idim=0; idim<Ndim;idim++){
         matrices[idim]=M.deepCopy();
-        for (int jdim=0; jdim<_Ndim+1;jdim++)
+        for (int jdim=0; jdim<Ndim+1;jdim++)
                        matrices[idim](jdim,idim) = values[jdim] ;
     }
 
-       Vector result(_Ndim);
-    for (int idim=0; idim<_Ndim;idim++)
+       Vector result(Ndim);
+    for (int idim=0; idim<Ndim;idim++)
         result[idim] = matrices[idim].determinant();
 
        return result;    
@@ -332,7 +323,7 @@ double StationaryDiffusionEquation::computeDiffusionMatrix(bool & stop)
 double StationaryDiffusionEquation::computeDiffusionMatrixFE(bool & stop){
        Cell Cj;
        string nameOfGroup;
-       double dn, coeff;
+       double coeff;
        MatZeroEntries(_A);
        VecZeroEntries(_b);
     
@@ -660,22 +651,13 @@ bool StationaryDiffusionEquation::iterateNewtonStep(bool &converged)
         VecAssemblyBegin(_deltaT);
         VecAssemblyEnd(  _deltaT);
 
-        if(!_FECalculation)
-            for(int i=0; i<_Nmailles; i++)
-            {
-                VecGetValues(_deltaT, 1, &i, &dTi);
-                VecGetValues(_Tk, 1, &i, &Ti);
-                if(_erreur_rel < fabs(dTi/Ti))
-                    _erreur_rel = fabs(dTi/Ti);
-            }
-        else
-            for(int i=0; i<_NunknownNodes; i++)
-            {
-                VecGetValues(_deltaT, 1, &i, &dTi);
-                VecGetValues(_Tk, 1, &i, &Ti);
-                if(_erreur_rel < fabs(dTi/Ti))
-                    _erreur_rel = fabs(dTi/Ti);
-            }
+               for(int i=0; i<_VV.getNumberOfElements(); i++)
+               {
+                       VecGetValues(_deltaT, 1, &i, &dTi);
+                       VecGetValues(_Tk, 1, &i, &Ti);
+                       if(_erreur_rel < fabs(dTi/Ti))
+                               _erreur_rel = fabs(dTi/Ti);
+               }
         stop=false;
         converged = (_erreur_rel <= _precision) ;//converged=convergence des iterations de Newton
     }
@@ -1005,10 +987,7 @@ StationaryDiffusionEquation::getEigenvectorsField(int nev, EPSWhich which, doubl
   MEDCoupling::DataArrayDouble * d = A.getEigenvectorsDataArrayDouble( nev, which, tol);
   Field my_eigenfield;
   
-  if(_FECalculation)
-    my_eigenfield = Field("Eigenvectors field", NODES, _mesh, nev);
-  else
-    my_eigenfield = Field("Eigenvectors field", CELLS, _mesh, nev);
+    my_eigenfield = Field("Eigenvectors field", _VV.getTypeOfField(), _mesh, nev);
 
   my_eigenfield.setFieldByDataArrayDouble(d);
   

diff --git a/CoreFlows/Models/src/TransportEquation.cxx b/CoreFlows/Models/src/TransportEquation.cxx
index 8db4423936b01f6184f6ef9f6c61000f813ffda4..cdd54132dcfb447ce8c6f3b66376b7e3313b46f1 100755 (executable)
--- a/CoreFlows/Models/src/TransportEquation.cxx
+++ b/CoreFlows/Models/src/TransportEquation.cxx
@@ -58,6 +58,8 @@ void TransportEquation::initialize()
 {
        if(!_initialDataSet)
                throw CdmathException("TransportEquation::initialize() set initial data first");
+       else if (_VV.getTypeOfField() != CELLS)
+               throw CdmathException("TransportEquation::initialize() Initial data should be a field on CELLS, not NODES, neither FACES");
        else
                cout<<"Initialising the transport of a fluid enthalpy"<<endl;
        /**************** Field creation *********************/