From: michael <michael@localhost.localdomain>
Date: Sat, 28 Nov 2020 23:08:34 +0000 (+0100)
Subject: Working on linear elasticity
X-Git-Tag: V9_7_0~129
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=5e1778c4d8442850f5c518867060186780455172;p=tools%2Fsolverlab.git

Working on linear elasticity
---

diff --git a/CoreFlows/Models/inc/LinearElasticityModel.hxx b/CoreFlows/Models/inc/LinearElasticityModel.hxx
index 741a306..81aeed2 100755
--- a/CoreFlows/Models/inc/LinearElasticityModel.hxx
+++ b/CoreFlows/Models/inc/LinearElasticityModel.hxx
@@ -59,14 +59,12 @@ public :
 
 	void setConstantDensity(double rho) { _rho=rho; }
 	void setDensityField(Field densityField) { _densityField=densityField; _densityFieldSet=true;}
-	void setLameCoefficient(double lambda, double mu) { _lambda = lambda; _mu = mu;}
+	void setLameCoefficients(double lambda, double mu) { _lambda = lambda; _mu = mu;}
 	void setYoungAndPoissonModuli(double E, double nu) { _lambda = E*nu/(1+nu)/(1-2*nu); _mu = E/2/(1+nu);}
 	void setGravity(Vector gravite ) { _gravity=gravite; }
 
     void setMesh(const Mesh &M);
-    void setFileName(string fileName){
-	_fileName = fileName;
-    }
+    void setFileName(string fileName){	_fileName = fileName;    }
     bool solveStationaryProblem();
     Field getOutputDisplacementField();
     
diff --git a/CoreFlows/Models/src/LinearElasticityModel.cxx b/CoreFlows/Models/src/LinearElasticityModel.cxx
index a81eebf..bdcd35c 100755
--- a/CoreFlows/Models/src/LinearElasticityModel.cxx
+++ b/CoreFlows/Models/src/LinearElasticityModel.cxx
@@ -221,11 +221,10 @@ void LinearElasticityModel::initialize()
 	PCSetType(_pc, _pctype);
 
     //Checking whether all boundary conditions are Neumann boundary condition
-    //if(_FECalculation) _onlyNeumannBC = _NdirichletNodes==0;
     if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
     {
-        map<string, LimitFieldStationaryDiffusion>::iterator it = _limitField.begin();
-        while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion)
+        map<string, LimitFieldLinearElasticity>::iterator it = _limitField.begin();
+        while(it != _limitField.end() and (it->second).bcType == NeumannLinearElasticity)
             it++;
         _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);//what if _limitField.size()==0 ???
     }
@@ -235,11 +234,6 @@ void LinearElasticityModel::initialize()
         else
             _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
 
-    //Checking whether all boundaries are Neumann boundaries
-    map<string, LimitFieldLinearElasticity>::iterator it = _limitField.begin();
-    while(it != _limitField.end() and (it->second).bcType == NeumannLinearElasticity)
-        it++;
-    _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);
     //If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
     if(_onlyNeumannBC)
     {
@@ -389,9 +383,9 @@ double LinearElasticityModel::computeStiffnessMatrixFE(bool & stop){
                 {
                     j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
                     if( _neumannValuesSet )
-                        coeff =Fi.getMeasure()/_Ndim*_neumannBoundaryValues[Fi.getNodeId(j)];
+                        coeff =Fi.getMeasure()/(_Ndim+1)*_neumannBoundaryValues[Fi.getNodeId(j)];
                     else    
-                        coeff =Fi.getMeasure()/_Ndim*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalForce;
+                        coeff =Fi.getMeasure()/(_Ndim+1)*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalForce;
                     VecSetValue(_b, j_int, coeff, ADD_VALUES);
                 }
             }
@@ -414,7 +408,7 @@ double LinearElasticityModel::computeStiffnessMatrixFV(bool & stop){
 	Cell Cell1,Cell2;
 	string nameOfGroup;
 	double inv_dxi, inv_dxj;
-	double barycenterDistance;
+	double barycenterDistance,coeff;
 	Vector normale(_Ndim);
 	double dn;
 	PetscInt idm, idn;
@@ -464,7 +458,12 @@ double LinearElasticityModel::computeStiffnessMatrixFV(bool & stop){
             {
                 nameOfGroup = Fj.getGroupName();
     
-                if (_limitField[nameOfGroup].bcType==NeumannLinearElasticity){//Nothing to do
+                if (_limitField[nameOfGroup].bcType==NeumannLinearElasticity){
+                    if( _neumannValuesSet )
+                        coeff =Fi.getMeasure()*_neumannBoundaryValues[j];
+                    else    
+                        coeff =Fi.getMeasure()*_limitField[nameOfGroup].normalForce;
+                    VecSetValue(_b,idm,    coeff, ADD_VALUES);
                 }
                 else if(_limitField[nameOfGroup].bcType==DirichletLinearElasticity){
                     barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());