Working on linear elasticity

diff --git a/CoreFlows/Models/inc/LinearElasticityModel.hxx b/CoreFlows/Models/inc/LinearElasticityModel.hxx
index 741a306e7c0e866caed2b7fb832369130e0eeaed..81aeed208cb6cd7d30781201508456f4c7a209bb 100755 (executable)
--- 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 a81eebfa4415b0b337faacab35a8aafd68625c05..bdcd35c313ed2ea208d54dc2fe6e2c2313ccf594 100755 (executable)
--- 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());