Commit 067d2102 authored by malossi's avatar malossi
Browse files

Added a method to get the value of the lagrange multiplier associated with a...

Added a method to get the value of the lagrange multiplier associated with a flow rate imposed on a boundary identified by a flag.

The LM can be used to get the correct value of the normal stress: (sigma dot n) dot n
parent 35266c73
......@@ -228,12 +228,29 @@ public:
// compute average pressure on a boundary face with given flag and a given solution
Real pressure(const EntityFlag& flag, const vector_type& sol);
//! Get the lagrange multiplier related to a flux imposed on a given part of the boundary.
/*!
* @param Flag flag of the boundary face associated with the flux and the Lagrange multiplier we want.
* @param BC BChandler containing the boundary conditions of the problem.
* @return Lagrange multiplier
*/
Real LagrangeMultiplier( const EntityFlag& Flag, bchandler_raw_type& BC );
// compute flux on a boundary face with given flag
Real flux(const EntityFlag& flag);
// compute average pressure on a boundary face with given flag
Real pressure(const EntityFlag& flag);
//! Get the lagrange multiplier related to a flux imposed on a given part of the boundary.
/*!
* @param Flag flag of the boundary face associated with the flux and the Lagrange multiplier we want.
* @param BC BChandler containing the boundary conditions of the problem.
* @param Solution vector containing the solution of the problem (and also the Lagrange multipliers at the end).
* @return Lagrange multiplier
*/
Real LagrangeMultiplier( const EntityFlag& Flag, bchandler_raw_type& BC, const vector_type& Solution );
// return the density and the viscosity of the fluid
Real density() const { return M_data.density(); }
Real viscosity() const { return M_data.viscosity(); }
......@@ -1318,7 +1335,6 @@ Oseen<Mesh, SolverType>::flux(const EntityFlag& flag)
return flux(flag, M_sol);
}
//! Computes the pressure on a given part of the boundary
template<typename Mesh, typename SolverType> Real
Oseen<Mesh, SolverType>::pressure(const EntityFlag& flag)
......@@ -1326,6 +1342,12 @@ Oseen<Mesh, SolverType>::pressure(const EntityFlag& flag)
return pressure(flag, M_sol);
}
template<typename Mesh, typename SolverType> Real
Oseen<Mesh, SolverType>::LagrangeMultiplier( const EntityFlag& Flag, bchandler_raw_type& BC )
{
return LagrangeMultiplier( Flag, BC, M_sol );
}
//! Computes the flux on a given part of the boundary
template<typename Mesh, typename SolverType> Real
Oseen<Mesh, SolverType>::flux(const EntityFlag& flag, const vector_type& sol)
......@@ -1349,6 +1371,23 @@ Oseen<Mesh, SolverType>::pressure(const EntityFlag& flag, const vector_type& sol
return M_post_proc->average(press, flag, 1)[0];
}
template<typename Mesh, typename SolverType> Real
Oseen<Mesh, SolverType>::LagrangeMultiplier( const EntityFlag& Flag, bchandler_raw_type& BC, const vector_type& Solution )
{
// Create a list of Flux BCName
std::vector< BCName > FluxBCVector = BC.getBCWithType( Flux );
BCName FluxBCName = BC.GetBCWithFlag( Flag ).name();
// Find the index associated to the correct Lagrange multiplier
for ( UInt lmIndex(0); lmIndex < static_cast <UInt> ( FluxBCVector.size() ); ++lmIndex )
if ( FluxBCName.compare( FluxBCVector[lmIndex] ) == 0 )
return Solution[3 * M_uFESpace.dof().numTotalDof() + M_pFESpace.dof().numTotalDof() + 1 + lmIndex];
// If lmIndex has not been found a warning message is printed
std::cout << "!!! Warning - Lagrange multiplier for Flux BC not found!" << std::endl;
return 0;
}
//! Computes the area on a given part of the boundary
template<typename Mesh, typename SolverType> Real
Oseen<Mesh, SolverType>::area(const EntityFlag& flag) {
......
......@@ -100,6 +100,20 @@ public:
Real GetLinearFlux ( const EntityFlag& flag );
Real GetLinearPressure( const EntityFlag& flag );
//! Get the lagrange multiplier related to a flux imposed on a given part of the boundary.
/*!
* @param Flag flag of the boundary face associated with the flux and the Lagrange multiplier we want.
* @param BC BChandler containing the boundary conditions of the problem.
* @return Lagrange multiplier
*/
Real LinearLagrangeMultiplier( const EntityFlag& Flag, bchandler_raw_type& BC );
//! Get the solution of the Shape Derivative problem.
/*!
* @return vector containing the solution of the Shape Derivative problem.
*/
const vector_type& LinearSolution() const { return M_linSol; }
void updateShapeDerivatives(
matrix_type& matrNoBC,
double& alpha,
......@@ -113,7 +127,7 @@ public:
bool convectiveTermDerivative=false);
void updateRhsLinNoBC( const vector_type& rhs){M_rhsLinNoBC=rhs;}
bool const stab(){return this->M_stab;}
bool stab() { return this->M_stab; }
private:
......@@ -254,18 +268,27 @@ OseenShapeDerivative<Mesh, SolverType>::
}
template<typename Mesh, typename SolverType> Real
template<typename Mesh, typename SolverType>
Real
OseenShapeDerivative<Mesh, SolverType>::GetLinearFlux( const EntityFlag& flag )
{
return flux( flag, M_linSol );
}
template<typename Mesh, typename SolverType> Real
template<typename Mesh, typename SolverType>
Real
OseenShapeDerivative<Mesh, SolverType>::GetLinearPressure( const EntityFlag& flag )
{
return pressure( flag, M_linSol );
}
template<typename Mesh, typename SolverType>
Real
OseenShapeDerivative<Mesh, SolverType>::LinearLagrangeMultiplier( const EntityFlag& Flag, bchandler_raw_type& BC )
{
return LagrangeMultiplier( Flag, BC, M_linSol );
}
template<typename Mesh, typename SolverType>
void OseenShapeDerivative<Mesh, SolverType>::setUp( const GetPot& dataFile )
{
......@@ -477,7 +500,7 @@ OseenShapeDerivative<Mesh, SolverType>::updateLinearSystem( const matrix_type& /
M_elvec_du.showMe(std::cout);
std::cout << "fin ====================" << std::endl;
*/
// assembling presssure right hand side
// assembling pressure right hand side
assembleVector( rhsLinNoBC, M_elvec_dp, this->M_pFESpace.fe(), this->M_pFESpace.dof(), 0, nbCompU*this->dim_u() );
// loop on velocity components
......@@ -522,7 +545,7 @@ OseenShapeDerivative<Mesh, SolverType>::updateShapeDerivatives( matrix_type& M_m
Chrono chrono;
int nbCompU = nDimensions;
UInt nbCompU = nDimensions;
// M_rhsLinNoBC = sourceVec;//which is usually zero
......@@ -584,7 +607,7 @@ OseenShapeDerivative<Mesh, SolverType>::updateShapeDerivatives( matrix_type& M_m
{
UInt iloc = this->M_uFESpace.fe().patternFirst( k ); // iloc = k
for ( int ic = 0; ic < nbCompU; ++ic )
for ( UInt ic = 0; ic < nbCompU; ++ic )
{
UInt ig = this->M_uFESpace.dof().localToGlobal( i, iloc + 1 ) + ic * this->dim_u();
......@@ -694,9 +717,7 @@ OseenShapeDerivative<Mesh, SolverType>::updateShapeDerivatives( matrix_type& M_m
chrono.stop();
this->M_Displayer.leaderPrintMax("done in ", chrono.diff() );
}
//#endif
}
#endif
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment