@@ -41,7 +41,7 @@ BEGIN_FECORE_CLASS(FEBiphasicSolidDomain, FESolidDomain)
4141 ADD_PARAMETER (m_secant_stress, " secant_stress" );
4242 ADD_PARAMETER (m_secant_tangent, " secant_tangent" );
4343 ADD_PARAMETER (m_secant_perm_tangent, " secant_permeability_tangent" );
44- ADD_PARAMETER (m_auto_pressure_stab, " auto_pressure_stabilization" );
44+ ADD_PARAMETER (m_auto_pressure_stab, " auto_pressure_stabilization" )->setEnums( " off \0 on \0 adaptive \0 " ) ;
4545END_FECORE_CLASS ();
4646
4747// -----------------------------------------------------------------------------
@@ -51,7 +51,7 @@ FEBiphasicSolidDomain::FEBiphasicSolidDomain(FEModel* pfem) : FESolidDomain(pfem
5151 m_secant_stress = false ;
5252 m_secant_tangent = false ;
5353 m_secant_perm_tangent = false ;
54- m_auto_pressure_stab = false ;
54+ m_auto_pressure_stab = 0 ;
5555
5656 if (pfem)
5757 {
@@ -138,6 +138,16 @@ void FEBiphasicSolidDomain::PreSolveUpdate(const FETimeInfo& timeInfo)
138138 mp.Update (timeInfo);
139139 }
140140 }
141+
142+ // see if we update the stabilization factor
143+ if (m_auto_pressure_stab > 1 )
144+ {
145+ if (!CalcAutoPressureStabilization ())
146+ {
147+ feLogError (" Failed to calculate automatic pressure stabilization factor." );
148+ return ;
149+ }
150+ }
141151}
142152
143153// -----------------------------------------------------------------------------
@@ -158,7 +168,7 @@ bool FEBiphasicSolidDomain::Init()
158168 // allocate nodal pressures
159169 m_nodePressure.resize (Nodes (), 0.0 );
160170
161- if (m_auto_pressure_stab)
171+ if (m_auto_pressure_stab > 0 )
162172 {
163173 if (!CalcAutoPressureStabilization ())
164174 {
@@ -1267,6 +1277,8 @@ void FEBiphasicSolidDomain::GetNodalPressures(vector<double>& data)
12671277
12681278bool FEBiphasicSolidDomain::CalcAutoPressureStabilization ()
12691279{
1280+ feLog (" Updating stabilization factor for biphasic solid domain '%s'\n " , GetName ().c_str ());
1281+
12701282 // make sure that the biphasic tau parameter is not zero
12711283 if (m_pMat->m_tau == 0.0 )
12721284 {
@@ -1282,7 +1294,9 @@ bool FEBiphasicSolidDomain::CalcAutoPressureStabilization()
12821294
12831295 FEMesh& mesh = *GetMesh ();
12841296
1285- for (int i = 0 ; i < Elements (); ++i)
1297+ int NE = Elements ();
1298+ #pragma omp parallel for shared(EL, EEL, mesh)
1299+ for (int i = 0 ; i < NE ; ++i)
12861300 {
12871301 FESolidElement& el = Element (i);
12881302
@@ -1304,7 +1318,6 @@ bool FEBiphasicSolidDomain::CalcAutoPressureStabilization()
13041318 {
13051319 FEMaterialPoint& mp = *el.GetMaterialPoint (n);
13061320 FEBiphasicMaterialPoint* bpt = (mp.ExtractData <FEBiphasicMaterialPoint>());
1307- if (bpt == nullptr ) return false ;
13081321
13091322 // get the tangent stiffness
13101323 tens4dmm C = m_pMat->Tangent (mp);
0 commit comments