From 5e1778c4d8442850f5c518867060186780455172 Mon Sep 17 00:00:00 2001 From: michael Date: Sun, 29 Nov 2020 00:08:34 +0100 Subject: [PATCH] Working on linear elasticity --- .../Models/inc/LinearElasticityModel.hxx | 6 ++--- .../Models/src/LinearElasticityModel.cxx | 23 +++++++++---------- 2 files changed, 13 insertions(+), 16 deletions(-) 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::iterator it = _limitField.begin(); - while(it != _limitField.end() and (it->second).bcType == NeumannStationaryDiffusion) + map::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::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()); -- 2.39.2