_NdirichletNodes=0;
_NunknownNodes=0;
_dirichletValuesSet=false;
+ _neumannValuesSet=false;
//Linear solver data
_precision=1.e-6;
_runLogFile->close();
throw CdmathException("Missing boundary value");
}
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoTypeSpecified)
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==NoneBCLinearElasticity)
{
cout<<"!!! No boundary condition set for boundary node " << _boundaryNodeIds[i]<< endl;
*_runLogFile<< "!!!No boundary condition set for boundary node " << _boundaryNodeIds[i]<<endl;
_runLogFile->close();
throw CdmathException("Missing boundary condition");
}
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==Dirichlet)
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType==DirichletLinearElasticity)
_dirichletNodeIds.push_back(_boundaryNodeIds[i]);
- else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=Neumann)
+ else if(_limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType!=NeumannLinearElasticity)
{
cout<<"!!! Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<< endl;
- cout<<"!!! Accepted boundary conditions are Dirichlet "<< Dirichlet <<" and Neumann "<< Neumann << endl;
+ cout<<"!!! Accepted boundary conditions are Dirichlet "<< DirichletLinearElasticity <<" and Neumann "<< NeumannLinearElasticity << endl;
*_runLogFile<< "Wrong boundary condition "<< _limitField[_mesh.getNode(_boundaryNodeIds[i]).getGroupName()].bcType<< " set for boundary node " << _boundaryNodeIds[i]<<endl;
- *_runLogFile<< "Accepted boundary conditions are Dirichlet "<< Dirichlet <<" and Neumann "<< Neumann <<endl;
+ *_runLogFile<< "Accepted boundary conditions are Dirichlet "<< DirichletLinearElasticity <<" and Neumann "<< NeumannLinearElasticity <<endl;
_runLogFile->close();
throw CdmathException("Wrong boundary condition");
}
else
MatCreateSeqAIJ(PETSC_COMM_SELF, _NunknownNodes*_nVar, _NunknownNodes*_nVar, (1+_neibMaxNbNodes), PETSC_NULL, &_A);
- VecCreate(PETSC_COMM_SELF, &_displacments);
+ VecCreate(PETSC_COMM_SELF, &_displacements);
VecDuplicate(_displacements, &_b);//RHS of the linear system
KSPGetPC(_ksp, &_pc);
PCSetType(_pc, _pctype);
- //Checking whether all boundaries are Neumann boundaries
- map<string, LimitField>::iterator it = _limitField.begin();
- while(it != _limitField.end() and (it->second).bcType == Neumann)
- it++;
- _onlyNeumannBC = (it == _limitField.end() && _limitField.size()>0);
+ //Checking whether all boundary conditions are Neumann boundary condition
+ if(!_neumannValuesSet)//Boundary conditions set via LimitField structure
+ {
+ 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 ???
+ }
+ else
+ if(_FECalculation)
+ _onlyNeumannBC = _neumannBoundaryValues.size()==_NboundaryNodes;
+ else
+ _onlyNeumannBC = _neumannBoundaryValues.size()==_mesh.getBoundaryFaceIds().size();
+
//If only Neumann BC, then matrix is singular and solution should be sought in space of mean zero vectors
if(_onlyNeumannBC)
{
return result;
}
-double LinearElasticityModel::computeDiffusionMatrix(bool & stop)
+double LinearElasticityModel::computeStiffnessMatrix(bool & stop)
{
double result;
if(_FECalculation)
- result=computeDiffusionMatrixFE(stop);
+ result=computeStiffnessMatrixFE(stop);
else
- result=computeDiffusionMatrixFV(stop);
+ result=computeStiffnessMatrixFV(stop);
if(_verbose or _system)
MatView(_A,PETSC_VIEWER_STDOUT_SELF);
return result;
}
-double LinearElasticityModel::computeDiffusionMatrixFE(bool & stop){
+double LinearElasticityModel::computeStiffnessMatrixFE(bool & stop){
Cell Cj;
string nameOfGroup;
- double dn;
+ double dn, coeff;
MatZeroEntries(_A);
VecZeroEntries(_b);
if( _dirichletValuesSet )//New way of storing BC
valuesBorder[kdim]=_dirichletBoundaryValues[it->second];
else //old way of storing BC
- valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].Displacements;
+ valuesBorder[kdim]=_limitField[_mesh.getNode(nodeIds[kdim]).getGroupName()].displacement;
}
else
valuesBorder[kdim]=Vector(_Ndim);
}
}
+ //Calcul de la contribution de la condition limite de Neumann au second membre
+ if( _NdirichletNodes !=_NboundaryNodes)
+ {
+ vector< int > boundaryFaces = _mesh.getBoundaryFaceIds();
+ int NboundaryFaces=boundaryFaces.size();
+ for(int i = 0; i< NboundaryFaces ; i++)//On parcourt les faces du bord
+ {
+ Face Fi = _mesh.getFace(i);
+ for(int j = 0 ; j<_Ndim ; j++)//On parcourt les noeuds de la face
+ {
+ if(find(_dirichletNodeIds.begin(),_dirichletNodeIds.end(),Fi.getNodeId(j))==_dirichletNodeIds.end())//node j is an Neumann BC node (not a Dirichlet BC node)
+ {
+ j_int=unknownNodeIndex(Fi.getNodeId(j), _dirichletNodeIds);//indice du noeud j en tant que noeud inconnu
+ if( _neumannValuesSet )
+ coeff =Fi.getMeasure()/(_Ndim+1)*_neumannBoundaryValues[Fi.getNodeId(j)];
+ else
+ coeff =Fi.getMeasure()/(_Ndim+1)*_limitField[_mesh.getNode(Fi.getNodeId(j)).getGroupName()].normalForce;
+ VecSetValue(_b, j_int, coeff, ADD_VALUES);
+ }
+ }
+ }
+ }
+
MatAssemblyBegin(_A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(_A, MAT_FINAL_ASSEMBLY);
VecAssemblyBegin(_b);
return INFINITY;
}
-double LinearElasticityModel::computeDiffusionMatrixFV(bool & stop){
+double LinearElasticityModel::computeStiffnessMatrixFV(bool & stop){
long nbFaces = _mesh.getNumberOfFaces();
Face Fj;
Cell Cell1,Cell2;
string nameOfGroup;
double inv_dxi, inv_dxj;
- double barycenterDistance;
+ double barycenterDistance,coeff;
Vector normale(_Ndim);
double dn;
PetscInt idm, idn;
{
nameOfGroup = Fj.getGroupName();
- if (_limitField[nameOfGroup].bcType==Neumann){//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==Dirichlet){
+ else if(_limitField[nameOfGroup].bcType==DirichletLinearElasticity){
barycenterDistance=Cell1.getBarryCenter().distance(Fj.getBarryCenter());
MatSetValue(_A,idm,idm,dn*inv_dxi/barycenterDistance , ADD_VALUES);
- VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].T, ADD_VALUES);
+ VecSetValue(_b,idm, dn*inv_dxi/barycenterDistance*_limitField[nameOfGroup].displacement, ADD_VALUES);
}
else {
stop=true ;
- cout<<"!!!!!!!!!!!!!!! Error LinearElasticityModel::computeDiffusionMatrixFV !!!!!!!!!!"<<endl;
+ cout<<"!!!!!!!!!!!!!!! Error LinearElasticityModel::computeStiffnessMatrixFV !!!!!!!!!!"<<endl;
cout<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
- cout<<"Accepted boundary conditions are Neumann "<<Neumann<< " and Dirichlet "<<Dirichlet<<endl;
+ cout<<"Accepted boundary conditions are Neumann "<<NeumannLinearElasticity<< " and Dirichlet "<<DirichletLinearElasticity<<endl;
*_runLogFile<<"!!!!!! Boundary condition not accepted for boundary named !!!!!!!!!!"<<nameOfGroup<< ", _limitField[nameOfGroup].bcType= "<<_limitField[nameOfGroup].bcType<<endl;
_runLogFile->close();
throw CdmathException("Boundary condition not accepted");
}
else
{
- *_runLogFile<<"LinearElasticityModel::computeDiffusionMatrixFV(): incompatible number of cells around a face"<<endl;
- throw CdmathException("LinearElasticityModel::computeDiffusionMatrixFV(): incompatible number of cells around a face");
+ *_runLogFile<<"LinearElasticityModel::computeStiffnessMatrixFV(): incompatible number of cells around a face"<<endl;
+ throw CdmathException("LinearElasticityModel::computeStiffnessMatrixFV(): incompatible number of cells around a face");
}
}
if(_conditionNumber)
KSPSetComputeEigenvalues(_ksp,PETSC_TRUE);
- KSPSolve(_ksp, _b, displacements);
+ KSPSolve(_ksp, _b, _displacements);
KSPConvergedReason reason;
KSPGetConvergedReason(_ksp,&reason);
*_runLogFile<< "Finite elements method"<< endl<<endl;
}
- computeDiffusionMatrix( stop);
+ computeStiffnessMatrix( stop);
if (stop){
cout << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
*_runLogFile << "Error : failed computing diffusion matrix, stopping calculation"<< endl;
_dirichletValuesSet=true;
_dirichletBoundaryValues=dirichletBoundaryValues;
}
+void
+LinearElasticityModel::setNeumannValues(map< int, double> neumannBoundaryValues)
+{
+ _neumannValuesSet=true;
+ _neumannBoundaryValues=neumannBoundaryValues;
+}
double
LinearElasticityModel::getConditionNumber(bool isSingular, double tol) const
