From 5aca970ac7b602bea93a1ec1d45632abca192a8f Mon Sep 17 00:00:00 2001 From: Wenqing Wang Date: Fri, 21 Dec 2018 14:25:28 +0100 Subject: [PATCH] [Excavation] Fixed a bug in the excavation modelling with HM processes --- FEM/fem_ele_vec.cpp | 169 +++++++++++++++----------------------------- FEM/rf_pcs.cpp | 21 ++++-- 2 files changed, 73 insertions(+), 117 deletions(-) diff --git a/FEM/fem_ele_vec.cpp b/FEM/fem_ele_vec.cpp index 5c56ab02e..e34d77c76 100644 --- a/FEM/fem_ele_vec.cpp +++ b/FEM/fem_ele_vec.cpp @@ -1026,7 +1026,7 @@ void CFiniteElementVec::LocalAssembly(const int update) } else #endif //#if !defined(USE_PETSC) // && !defined(other parallel libs)//03.3012. - //WW + // WW { #if defined(USE_PETSC) // || defined (other parallel solver lib). 04.2012 WW // TODO @@ -1056,7 +1056,7 @@ void CFiniteElementVec::LocalAssembly(const int update) } else #endif //#if !defined(USE_PETSC) // && !defined(other parallel libs)//03.3012. - //WW + // WW #if !defined(USE_PETSC) // && !defined(other parallel libs)//03.3012. WW { for (int i = 0; i < nnodesHQ; i++) @@ -1637,71 +1637,6 @@ void CFiniteElementVec::GlobalAssembly_RHS() if (dynamic) Residual = true; - bool onExBoundaryState[max_nnodes_QE_3D] = {false}; - if (excavation) - { - int valid = 0; - excavation = true; - bool onExBoundary = false; - CNode* node; - CElem* elem; - CSolidProperties* smat_e; - - for (int i = 0; i < nnodesHQ; i++) - { - node = MeshElement->nodes[i]; - onExBoundary = false; - const std::size_t n_elements(node->getConnectedElementIDs().size()); - for (std::size_t j = 0; j < n_elements; j++) - { - elem = - pcs->m_msh->ele_vector[node->getConnectedElementIDs()[j]]; - if (!elem->GetMark()) - continue; - - smat_e = msp_vector[elem->GetPatchIndex()]; - if (smat_e->excavation > 0) - { - if (fabs(GetCurveValue(smat_e->excavation, 0, aktuelle_zeit, - &valid) - - 1.0) < DBL_MIN) - { - onExBoundary = true; - break; - } - } - else if (pcs->ExcavMaterialGroup > -1) - { - double const* ele_center(elem->GetGravityCenter()); - const double expected_coordinate = - GetCurveValue(pcs->ExcavCurve, 0, aktuelle_zeit, - &valid) + - pcs->ExcavBeginCoordinate; - const double max_excavation_range = std::max( - expected_coordinate, pcs->ExcavBeginCoordinate); - if (ele_center[pcs->ExcavDirection] > max_excavation_range) - { - onExBoundary = true; - break; - } - else if (elem->GetPatchIndex() != - static_cast(pcs->ExcavMaterialGroup)) - { - onExBoundary = true; - break; - } - } - else - { - onExBoundary = true; - break; - } - } - if (onExBoundary) - onExBoundaryState[i] = 1; - } - } - if (Residual) { const double biot = smat->biot_const; @@ -1714,14 +1649,6 @@ void CFiniteElementVec::GlobalAssembly_RHS() case FiniteElement::LIQUID_FLOW: for (int i = 0; i < nnodes; i++) nodal_pore_p[i] = LoadFactor * _nodal_p1[i]; - if (excavation) - { - for (int i = 0; i < nnodes; i++) - { - if (onExBoundaryState[i] == 1) // WX:02.2013 - nodal_pore_p[i] = 0.0; - } - } if (pcs->Neglect_H_ini == 2) { for (int i = 0; i < nnodes; i++) @@ -1754,14 +1681,6 @@ void CFiniteElementVec::GlobalAssembly_RHS() nodal_pore_p[i] = LoadFactor * _nodal_S[i] * _nodal_p1[i]; } - if (excavation) - { - for (int i = 0; i < nnodes; i++) - { - if (onExBoundaryState[i] == 1) // WX:02.2013 - nodal_pore_p[i] = 0.0; - } - } if (pcs->Neglect_H_ini == 2) { for (int i = 0; i < nnodes; i++) @@ -1790,15 +1709,6 @@ void CFiniteElementVec::GlobalAssembly_RHS() smat->getBishopCoefficient(S_e, val_n) * val_n; } - if (excavation) - { - for (int i = 0; i < nnodes; i++) - { - if (onExBoundaryState[i] == 1) - nodal_pore_p[i] = 0.0; - } - } - if (pcs->Neglect_H_ini == 2) { for (int i = 0; i < nnodes; i++) @@ -1834,15 +1744,6 @@ void CFiniteElementVec::GlobalAssembly_RHS() nodal_pore_p[i] = pore_p * LoadFactor; } - if (excavation) - { - for (int i = 0; i < nnodes; i++) - { - if (onExBoundaryState[i] == 1) // WX:02.2013 - nodal_pore_p[i] = 0.; - } - } - if (pcs->Neglect_H_ini == 2) { for (int i = 0; i < nnodes; i++) @@ -1873,15 +1774,6 @@ void CFiniteElementVec::GlobalAssembly_RHS() nodal_pore_p[i] = pore_p * LoadFactor; } - if (excavation) - { - for (int i = 0; i < nnodes; i++) - { - if (onExBoundaryState[i] == 1) // WX:02.2013 - nodal_pore_p[i] = 0.; - } - } - if (pcs->Neglect_H_ini == 2) { for (int i = 0; i < nnodes; i++) @@ -1947,9 +1839,64 @@ void CFiniteElementVec::GlobalAssembly_RHS() if (excavation) { + int valid = 0; + excavation = true; + bool onExBoundary = false; + CNode* node; + CElem* elem; + CSolidProperties* smat_e; + for (int i = 0; i < nnodesHQ; i++) { - if (!onExBoundaryState[i]) + node = MeshElement->nodes[i]; + onExBoundary = false; + const std::size_t n_elements(node->getConnectedElementIDs().size()); + for (std::size_t j = 0; j < n_elements; j++) + { + elem = + pcs->m_msh->ele_vector[node->getConnectedElementIDs()[j]]; + if (!elem->GetMark()) + continue; + + smat_e = msp_vector[elem->GetPatchIndex()]; + if (smat_e->excavation > 0) + { + if (fabs(GetCurveValue(smat_e->excavation, 0, aktuelle_zeit, + &valid) - + 1.0) < DBL_MIN) + { + onExBoundary = true; + break; + } + } + else if (pcs->ExcavMaterialGroup > -1) + { + double const* ele_center(elem->GetGravityCenter()); + const double expected_coordinate = + GetCurveValue(pcs->ExcavCurve, 0, aktuelle_zeit, + &valid) + + pcs->ExcavBeginCoordinate; + const double max_excavation_range = std::max( + expected_coordinate, pcs->ExcavBeginCoordinate); + if (ele_center[pcs->ExcavDirection] > max_excavation_range) + { + onExBoundary = true; + break; + } + else if (elem->GetPatchIndex() != + static_cast(pcs->ExcavMaterialGroup)) + { + onExBoundary = true; + break; + } + } + else + { + onExBoundary = true; + break; + } + } + if (!onExBoundary) { for (std::size_t j = 0; j < dim; j++) b_rhs[eqs_number[i] + NodeShift[j]] = 0.0; diff --git a/FEM/rf_pcs.cpp b/FEM/rf_pcs.cpp index d554cff78..0063cfa0c 100644 --- a/FEM/rf_pcs.cpp +++ b/FEM/rf_pcs.cpp @@ -19,9 +19,9 @@ solver and parellelisation of them 02/2008 PCH OpenMP parallelization for Lis matrix solver **************************************************************************/ -#include "FEMEnums.h" -#include "Output.h" -#include "MathTools.h" +#include "rf_pcs.h" + +#include // std::numeric_limits /*--------------------- MPI Parallel -------------------*/ #if defined(USE_MPI) || defined(USE_MPI_PARPROC) || defined(USE_MPI_REGSOIL) @@ -55,11 +55,14 @@ matrix solver // GEOLib #include "PointWithID.h" +#include "FEMEnums.h" +#include "Output.h" +#include "MathTools.h" + #include "PhysicalConstant.h" /* Objects */ #include "pcs_dm.h" -#include "rf_pcs.h" #include "rf_st_new.h" // ST //#include "rf_bc_new.h" // ST //#include "rf_mmp_new.h" // MAT @@ -4617,8 +4620,14 @@ bool CRFProcess::isPointInExcavatedDomain(double const* point, point[ExcavDirection]; max_excavation_range = std::max(expected_coordinate, ExcavBeginCoordinate); min_excavation_range = std::min(expected_coordinate, ExcavBeginCoordinate); - if (element_center_x_in_excavation_direction > min_excavation_range && - element_center_x_in_excavation_direction < max_excavation_range) + if ((element_center_x_in_excavation_direction > min_excavation_range || + (std::fabs(element_center_x_in_excavation_direction - + min_excavation_range) < + std::numeric_limits::epsilon())) &&( + (element_center_x_in_excavation_direction < max_excavation_range) || + (std::fabs(element_center_x_in_excavation_direction - + max_excavation_range) < + std::numeric_limits::epsilon()))) { return true; }