From 75cb8badcc5ed8dd26e8d0e2f1029f6e45fe9144 Mon Sep 17 00:00:00 2001
From: michael <michael@localhost.localdomain>
Date: Sat, 30 Oct 2021 23:01:54 +0200
Subject: [PATCH] Better handling of field support (CELLS vs NODES)

---
 .../inc/StationaryDiffusionEquation.hxx       |  2 +-
 CoreFlows/Models/src/ProblemCoreFlows.cxx     | 11 ++-
 CoreFlows/Models/src/ProblemFluid.cxx         |  9 ++-
 .../src/StationaryDiffusionEquation.cxx       | 73 +++++++------------
 CoreFlows/Models/src/TransportEquation.cxx    |  2 +
 5 files changed, 44 insertions(+), 53 deletions(-)

diff --git a/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx b/CoreFlows/Models/inc/StationaryDiffusionEquation.hxx
index 7810a39..7bd065e 100755
--- 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 c19cc8d..aea9a50 100755
--- 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 8f497d3..731b5a7 100755
--- 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 a4661a3..44ae333 100755
--- 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 8db4423..cdd5413 100755
--- 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 *********************/
-- 
2.39.2