diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index e58b40478..85da4dd30 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -207,6 +207,8 @@ set(CSRCS ustruct.h ustruct.cpp vtk_xml.h vtk_xml.cpp vtk_xml_parser.h vtk_xml_parser.cpp + ris.h ris.cpp + uris.h uris.cpp CepMod.h CepMod.cpp CepModAp.h CepModAp.cpp diff --git a/Code/Source/solver/CmMod.cpp b/Code/Source/solver/CmMod.cpp index a0d7cd620..d5a97b1a3 100644 --- a/Code/Source/solver/CmMod.cpp +++ b/Code/Source/solver/CmMod.cpp @@ -138,6 +138,12 @@ void cmType::bcast(const CmMod& cm_mod, Array& data, const std::string& MPI_Bcast(data.data(), data.size(), cm_mod::mpreal, cm_mod.master, com()); } +/// @brief bcast int array +void cmType::bcast(const CmMod& cm_mod, Array& data, const std::string& name) const +{ + MPI_Bcast(data.data(), data.size(), cm_mod::mpint, cm_mod.master, com()); +} + /// @brief bcast double Vector void cmType::bcast(const CmMod& cm_mod, Vector& data, const std::string& name) const { @@ -155,12 +161,6 @@ void cmType::bcast(const CmMod& cm_mod, int* data) const /// @brief bcast int Vector void cmType::bcast(const CmMod& cm_mod, Vector& data) const -{ - MPI_Bcast(data.data(), data.size(), cm_mod::mpint, cm_mod.master, com()); -} - -/// @brief bcast int array -void cmType::bcast(const CmMod& cm_mod, Array& data, const std::string& name) const { MPI_Bcast(data.data(), data.size(), cm_mod::mpint, cm_mod.master, com()); } \ No newline at end of file diff --git a/Code/Source/solver/ComMod.cpp b/Code/Source/solver/ComMod.cpp index d461b123a..692616534 100644 --- a/Code/Source/solver/ComMod.cpp +++ b/Code/Source/solver/ComMod.cpp @@ -61,6 +61,7 @@ ComMod::ComMod() pstEq = false; sstEq = false; ibFlag = false; + risFlag = false; } diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 3980cdbd1..da2ff5f90 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -161,6 +161,9 @@ class bcType // Robin: apply only in normal direction bool rbnN = false; + // Strong/Weak application of Dirichlet BC + int clsFlgRis = 0; + // Pre/Res/Flat/Para... boundary types // // This stores differnt BCs as bitwise values. @@ -194,6 +197,9 @@ class bcType // Robin: damping double c = 0.0; + // RIS0D: resistance + double resistance = 0.0; + // Penalty parameters for weakly applied Dir BC Vector tauB{0.0, 0.0}; //double tauB[2]; @@ -961,6 +967,15 @@ class mshType /// @brief IB: Mesh size parameter double dx = 0.0; + /// @brief RIS resistance value + double res = 0.0; + + /// @brief RIS projection tolerance + double tol = 0.0; + + /// @brief The volume of this mesh + double v = 0.0; + /// @breif ordering: node ordering for boundaries std::vector> ordering; @@ -1061,6 +1076,13 @@ class mshType /// @brief IB: tracers traceType trc; + /// @brief RIS: flags of whether elemets are adjacent to RIS projections + // std::vector eRIS; + Vector eRIS; + + /// @brief RIS: processor ids to change element partitions to + Vector partRIS; + /// @brief TET4 quadrature modifier double qmTET4 = (5.0+3.0*sqrt(5.0))/20.0; @@ -1115,6 +1137,9 @@ class eqType /// @brief IB: Number of possible outputs int nOutIB = 0; + /// @brief URIS: Number of possible outputs + int nOutURIS = 0; + /// @brief Number of domains int nDmn = 0; @@ -1203,6 +1228,9 @@ class eqType /// @brief IB: Outputs std::vector outIB; + /// @brief URIS: Outputs + std::vector outURIS; + /// @brief Body force associated with this equation std::vector bf; }; @@ -1400,6 +1428,113 @@ class ibType ibCommType cm; }; +/// @brief Data type for Resistive Immersed Surface +// +class risFaceType +{ + public: + + /// @brief Number of RIS surface + int nbrRIS = 0; + + /// @brief Count time steps where no check is needed + Vector nbrIter; + + /// @brief List of meshes, and faces connected. The first face is the + // proximal pressure's face, while the second is the distal one + Array3 lst; + + /// @brief Resistance value + Vector Res; + + /// @brief Flag closed surface active, the valve is considerd open initially + std::vector clsFlg; + + /// @brief Mean distal and proximal pressure (1: distal, 2: proximal) + Array meanP; + + /// @brief Mean flux on the RIS surface + Vector meanFl; + + /// @brief Status RIS interface + std::vector status; +}; + +/// @brief Unfitted Resistive Immersed surface data type +// +class urisType +{ + public: + + // Name of the URIS instance. + std::string name; + + // Whether any file has been saved. + bool savedOnce = false; + + // Total number of IB nodes. + int tnNo = 0; + + // Number of IB meshes. + int nFa = 0; + + // Position coordinates (2D array: rows x columns). + Array x; + + // Displacement (new) (2D array). + Array Yd; + + // Default signed distance value away from the valve. + double sdf_default = 10.0; + + // Default distance value of the valve boundary (valve thickness). + double sdf_deps = 0.04; + + // Default distance value of the valve boundary when the valve is closed. + double sdf_deps_close = 0.25; + + // Displacements of the valve when it opens (3D array). + Array3 DxOpen; + + // Displacements of the valve when it closes (3D array). + Array3 DxClose; + + // Normal vector pointing in the positive flow direction (1D array). + Vector nrm; + + // Close flag. + bool clsFlg = true; + + // Iteration count. + int cnt = 1000000; + + // URIS: signed distance function of each node to the uris (1D array). + Vector sdf; + + // Mesh scale factor. + double scF = 1.0; + + // Mean pressure upstream. + double meanPU = 0.0; + + // Mean pressure downstream. + double meanPD = 0.0; + + // Relaxation factor to compute weighted averages of pressure values. + double relax_factor = 0.5; + + // Array to store the fluid mesh elements that the uris node is in (2D array). + Array elemId; + + // Array to count how many times a uris node is found in the fluid mesh of a processor (1D array). + Vector elemCounter; + + // Derived type variables + // IB meshes + std::vector msh; + +}; + /// @brief The ComMod class duplicates the data structures in the Fortran COMMOD module /// defined in MOD.f. /// @@ -1470,9 +1605,29 @@ class ComMod { /// @brief Postprocess step - convert bin to vtk bool bin2VTK = false; + /// @brief Whether any RIS surface is considered + bool risFlag = false; + + /// @brief Whether any one-sided RIS surface with 0D coupling is considered + bool ris0DFlag = false; + + /// @brief Whether any URIS surface is considered + bool urisFlag = false; + + /// @brief Whether the URIS surface is active + bool urisActFlag = false; + + /// @brief Number of URIS surfaces (uninitialized, to be set later) + int nUris; + + /// @brief URIS resistance + double urisRes; + + /// @brief URIS resistance when the valve is closed + double urisResClose; + /// @brief Whether to use precomputed state-variable solutions bool usePrecomp = false; - //----- int members -----// /// @brief Current domain @@ -1532,7 +1687,7 @@ class ComMod { /// @brief Total number of degrees of freedom per node int tDof = 0; - /// @brief Total number of nodes (total number of nodes on current processor across + /// @brief Total number of nodes (number of nodes on current proc across /// all meshes) int tnNo = 0; @@ -1542,6 +1697,9 @@ class ComMod { /// @brief Number of stress values to be stored int nsymd = 0; + /// @brief Nbr of iterations + int RisnbrIter = 0; + //----- double members -----// @@ -1601,6 +1759,18 @@ class ComMod { /// @brief IB: iblank used for immersed boundaries (1 => solid, 0 => fluid) Vector iblank; + /// @brief TODO: for now, better to organize these within a class + struct Array2D { + // std::vector> map; + Array map; + }; + + /// @brief RIS mapping array, with local (mesh) enumeration + std::vector risMapList; + + /// @brief RIS mapping array, with global (total) enumeration + std::vector grisMapList; + /// @brief Old time derivative of variables (acceleration); known result at current time step Array Ao; @@ -1692,6 +1862,12 @@ class ComMod { /// @brief IB: Immersed boundary data structure ibType ib; + /// @brief risFace object + risFaceType ris; + + /// @brief unfitted RIS object + std::vector uris; + bool debug_active = false; Timer timer; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index ae55e6ba4..a818ab9c6 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -44,6 +44,8 @@ // - MeshParameters - Mesh parameter section // - EquationParameters - Equation parameter section // - ProjectionParameters - Projection parameter section +// - RISProjectionParameters - RIS Projection parameter section +// - URISMeshParameters - Mesh parameter section // // These section objects may also contain objects representing the sub-sections // defined for each section. @@ -174,6 +176,12 @@ void Parameters::read_xml(std::string file_name) // Set Add_equation values. set_equation_values(root_element); + + // Set RIS projection values. + set_RIS_projection_values(root_element); + + // Set RIS projection values. + set_URIS_mesh_values(root_element); } void Parameters::set_contact_values(tinyxml2::XMLElement* root_element) @@ -248,6 +256,40 @@ void Parameters::set_projection_values(tinyxml2::XMLElement* root_element) } } +void Parameters::set_RIS_projection_values(tinyxml2::XMLElement* root_element) +{ + auto add_RIS_proj_item = root_element->FirstChildElement(RISProjectionParameters::xml_element_name_.c_str()); + + while (add_RIS_proj_item) { + const char* RIS_proj_name; + auto result = add_RIS_proj_item->QueryStringAttribute("name", &RIS_proj_name); + + RISProjectionParameters* RIS_proj_params = new RISProjectionParameters(); + RIS_proj_params->name.set(std::string(RIS_proj_name)); + RIS_proj_params->set_values(add_RIS_proj_item); + RIS_projection_parameters.push_back(RIS_proj_params); + + add_RIS_proj_item = add_RIS_proj_item->NextSiblingElement(RISProjectionParameters::xml_element_name_.c_str()); + } +} + +void Parameters::set_URIS_mesh_values(tinyxml2::XMLElement* root_element) +{ + auto add_URIS_mesh_item = root_element->FirstChildElement(URISMeshParameters::xml_element_name_.c_str()); + + while (add_URIS_mesh_item) { + const char* URIS_mesh_name; + auto result = add_URIS_mesh_item->QueryStringAttribute("name", &URIS_mesh_name); + + URISMeshParameters* URIS_mesh_params = new URISMeshParameters(); + URIS_mesh_params->name.set(std::string(URIS_mesh_name)); + URIS_mesh_params->set_values(add_URIS_mesh_item); + URIS_mesh_parameters.push_back(URIS_mesh_params); + + add_URIS_mesh_item = add_URIS_mesh_item->NextSiblingElement(URISMeshParameters::xml_element_name_.c_str()); + } +} + ////////////////////////////////////////////////////////// // BodyForceParameters // ////////////////////////////////////////////////////////// @@ -415,6 +457,8 @@ BoundaryConditionParameters::BoundaryConditionParameters() set_parameter("Weakly_applied", false, !required, weakly_applied); set_parameter("Zero_out_perimeter", false, !required, zero_out_perimeter); + + set_parameter("Resistance", 1.e5, !required, resistance); } void BoundaryConditionParameters::print_parameters() @@ -2539,6 +2583,189 @@ void ProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) xml_util_set_parameters(ftpr, xml_elem, error_msg); } +////////////////////////////////////////////////////////// +// RISProjectionParameters // +////////////////////////////////////////////////////////// + +/// @brief Define the XML element name for mesh parameters. +const std::string RISProjectionParameters::xml_element_name_ = "Add_RIS_projection"; + +RISProjectionParameters::RISProjectionParameters() +{ + // A parameter that must be defined. + bool required = true; + + name = Parameter("name", "", required); + + set_parameter("Project_from_face", "", required, project_from_face); + set_parameter("Resistance", 1.e6, !required, resistance); + set_parameter("Projection_tolerance", 0.0, !required, projection_tolerance); +} + +void RISProjectionParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + using namespace tinyxml2; + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + + // Get the 'type' from the element. + const char* sname; + auto result = xml_elem->QueryStringAttribute("name", &sname); + if (sname == nullptr) { + throw std::runtime_error("No TYPE given in the XML element."); + } + name.set(std::string(sname)); + + using std::placeholders::_1; + using std::placeholders::_2; + + std::function ftpr = + std::bind( &RISProjectionParameters::set_parameter_value, *this, _1, _2); + + xml_util_set_parameters(ftpr, xml_elem, error_msg); +} + + +////////////////////////////////////////////////////////// +// URIS Mesh Parameters // +////////////////////////////////////////////////////////// +// [HZ] implemente URIS parameters here + +// Process parameters for the 'Add_URIS_mesh' XML element used for defining URIS mesh elements. + +/// @brief Define the XML element name for mesh parameters. +const std::string URISMeshParameters::xml_element_name_ = "Add_URIS_mesh"; + +URISMeshParameters::URISMeshParameters() +{ + bool required = true; + + // Mesh name from Add_mesh element. + name = Parameter("name", "", required); + + // Parameters under Add_mesh element. + // + set_parameter("Mesh_scale_factor", 1.0, !required, mesh_scale_factor); + set_parameter("Thickness", 0.04, !required, thickness); + set_parameter("Closed_thickness", 0.25, !required, close_thickness); + set_parameter("Resistance", 1.0e5, !required, resistance); + set_parameter("Closed_resistance", 1.0e5, !required, resistance_close); + set_parameter("Valve_starts_as_closed", true, !required, valve_starts_as_closed); + set_parameter("Positive_flow_normal_file_path", "", !required, positive_flow_normal_file_path); +} + +void URISMeshParameters::print_parameters() +{ + std::cout << std::endl; + std::cout << "---------------" << std::endl; + std::cout << "URIS Mesh Parameters" << std::endl; + std::cout << "---------------" << std::endl; + std::cout << name.name() << ": " << name.value() << std::endl; + + auto params_name_value = get_parameter_list(); + for (auto& [ key, value ] : params_name_value) { + std::cout << key << ": " << value << std::endl; + } + + for (auto& face : URIS_face_parameters) { + face->print_parameters(); + } +} + +void URISMeshParameters::set_values(tinyxml2::XMLElement* mesh_elem) +{ + using namespace tinyxml2; + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + auto item = mesh_elem->FirstChildElement(); + + while (item != nullptr) { + auto name = std::string(item->Value()); + + // Add_face sub-element. + if (name == URISFaceParameters::xml_element_name_) { + auto URIS_face_params = new URISFaceParameters(); + URIS_face_params->set_values(item); + URIS_face_parameters.push_back(URIS_face_params); + } else if (item->GetText() != nullptr) { + auto value = item->GetText(); + try { + set_parameter_value(name, value); + } catch (const std::bad_function_call& exception) { + throw std::runtime_error(error_msg + name + "'."); + } + } else { + throw std::runtime_error(error_msg + name + "'."); + } + + item = item->NextSiblingElement(); + } +} + + +////////////////////////////////////////////////////////// +// URIS Face Parameters // +////////////////////////////////////////////////////////// + +/// @brief Process parameters for the 'Add_URIS_face' XML element. +/// +/// Define the XML element name for face parameters. +const std::string URISFaceParameters::xml_element_name_ = "Add_URIS_face"; + +URISFaceParameters::URISFaceParameters() +{ + set_xml_element_name(xml_element_name_); + + // A parameter that must be defined. + bool required = true; + + name = Parameter("name", "", required); + + set_parameter("Face_file_path", "", !required, face_file_path); + set_parameter("Open_motion_file_path", "", !required, open_motion_file_path); + set_parameter("Close_motion_file_path", "", !required, close_motion_file_path); + + // set_parameter("End_nodes_face_file_path", "", !required, end_nodes_face_file_path); + // set_parameter("Quadrature_modifier_TRI3", (2.0/3.0), !required, quadrature_modifier_TRI3); +} + +void URISFaceParameters::print_parameters() +{ + std::cout << std::endl; + std::cout << "---------------" << std::endl; + std::cout << "URIS Face Parameters" << std::endl; + std::cout << "---------------" << std::endl; + std::cout << name.name() << ": " << name.value() << std::endl; + std::cout << face_file_path.name() << ": " << face_file_path.value() << std::endl; +} + +void URISFaceParameters::set_values(tinyxml2::XMLElement* face_elem) +{ + using namespace tinyxml2; + + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + const char* face_name; + auto result = face_elem->QueryStringAttribute("name", &face_name); + name.set(std::string(face_name)); + auto item = face_elem->FirstChildElement(); + + while (item != nullptr) { + auto name = std::string(item->Value()); + auto value = item->GetText(); + + if (value == nullptr) { + throw std::runtime_error(error_msg + name + "'."); + } + + try { + set_parameter_value(name, value); + } catch (const std::bad_function_call& exception) { + throw std::runtime_error(error_msg + name + "'."); + } + + item = item->NextSiblingElement(); + } +} + + ////////////////////////////////////////////////////////// // LinearAlgebraParameters // ////////////////////////////////////////////////////////// diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 5d3ebd410..b70b72246 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -79,6 +79,7 @@ template /// 2) Mesh (MeshParameters) /// 3) Equation (EquationParameters) /// 4) Projection (ProjectionParameters) +/// 5) RIS Projection (RIS ProjectionParameters) /// /// Each object contains methods to parse the XML file for the parameters defined for it. /// These section objects may also contain objects representing the sub-sections defined @@ -818,6 +819,8 @@ class BoundaryConditionParameters : public ParameterLists Parameter value; Parameter weakly_applied; Parameter zero_out_perimeter; + + Parameter resistance; }; /// @brief The OutputParameters class stores parameters for the @@ -1528,6 +1531,100 @@ class MeshParameters : public ParameterLists Parameter quadrature_modifier_TET4; }; +////////////////////////////////////////////////////////// +// Resistive immersed surfaces method // +////////////////////////////////////////////////////////// + +/// @brief The RISProjectionParameters class stores parameters for the +/// 'Add_RIS_projection' XML element used for RIS valve simulations. +/// \code {.xml} +/// +/// right_ris +/// 1.e6 +/// +/// \endcode +class RISProjectionParameters : public ParameterLists +{ + public: + RISProjectionParameters(); + + void set_values(tinyxml2::XMLElement* xml_elem); + + static const std::string xml_element_name_; + + Parameter name; + + Parameter project_from_face; + Parameter resistance; + Parameter projection_tolerance; +}; + +/// @brief The URISFaceParameters class is used to store parameters for the +/// 'Add_URIS_face' XML element. +class URISFaceParameters : public ParameterLists +{ + public: + URISFaceParameters(); + + void print_parameters(); + void set_values(tinyxml2::XMLElement* xml_elem); + + static const std::string xml_element_name_; + + Parameter name; // Name of the valve surface + + Parameter face_file_path; // File path for the valve surface + Parameter open_motion_file_path; // File path for the open motion of the valve + Parameter close_motion_file_path; // File path for the close motion of the valve + +}; + +/// @brief The URISMeshParameters class is used to store paramaters for the +/// 'Add_URIS_mesh' XML element. +/// +/// \code {.xml} +/// +/// +/// meshes/uris_face.vtu +/// meshes/uris_facemotion_open.dat +/// meshes/uris_facemotion_close.dat +/// +/// 1.0 +/// 0.25 +/// 1.0e5 +/// meshes/normal.dat +/// +/// \endcode +class URISMeshParameters : public ParameterLists +{ + public: + URISMeshParameters(); + + static const std::string xml_element_name_; + + void print_parameters(); + void set_values(tinyxml2::XMLElement* mesh_elem); + std::string get_name() const { return name.value(); }; + // std::string get_path() const { return mesh_file_path.value(); }; + + std::vector URIS_face_parameters; + + // Add_mesh name + Parameter name; // Name of the valve mesh + + // Parameters under Add_URIS_mesh + Parameter mesh_scale_factor; // Scale factor for the mesh + Parameter thickness; // Thickness of the valve + Parameter close_thickness; // Thickness of the valve when it is closed + Parameter resistance; // Resistance of the valve + Parameter resistance_close; // Resistance of the valve when it is closed + Parameter valve_starts_as_closed; // Whether the valve starts as closed + Parameter positive_flow_normal_file_path; // File path for the positive flow normal + +}; + + + /// @brief The Parameters class stores parameter values read in from a solver input file. class Parameters { @@ -1549,6 +1646,9 @@ class Parameters { void set_projection_values(tinyxml2::XMLElement* root_element); void set_svzerodsolver_interface_values(tinyxml2::XMLElement* root_element); + void set_RIS_projection_values(tinyxml2::XMLElement* root_element); + void set_URIS_mesh_values(tinyxml2::XMLElement* root_element); + // Objects representing each parameter section of XML file. ContactParameters contact_parameters; GeneralSimulationParameters general_simulation_parameters; @@ -1556,6 +1656,10 @@ class Parameters { std::vector equation_parameters; std::vector projection_parameters; PrecomputedSolutionParameters precomputed_solution_parameters; + + std::vector RIS_projection_parameters; + std::vector URIS_mesh_parameters; + }; #endif diff --git a/Code/Source/solver/all_fun.cpp b/Code/Source/solver/all_fun.cpp index 68c86ceb4..0a0f6dff2 100644 --- a/Code/Source/solver/all_fun.cpp +++ b/Code/Source/solver/all_fun.cpp @@ -295,6 +295,153 @@ global(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, const Arra return result; } +/////////////////////////////////////////////////////////////////////////// +/// Convert the vIntegM function from svFSI to c++. +/// @brief This routine integrated a scalar field over a particular domain. +/// +/// @param iM mesh index +/// @param s an array containing a scalar value for each node in the mesh +/// Replicates 'FUNCTION vIntegM(dId, s, l, u, pFlag)' defined in ALLFUN.f. +// +double integ(const ComMod& com_mod, const CmMod& cm_mod, int iM, const Array& s) +{ + using namespace consts; + + #define n_debug_integ_v + #ifdef debug_integ_v + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "vIntegM " << " "; + dmsg << "iM: " << iM; + #endif + + int nNo = s.ncols(); + int tnNo = com_mod.tnNo; + bool ibFlag = com_mod.ibFlag; + + if (nNo != tnNo) { + std::string msg = "Incompatible vector size in vIntegM in domain: "; + msg += std::to_string(iM); + throw std::runtime_error(msg); + } + + double result = 0.0; + int nsd = com_mod.nsd; + fsType fs; + + auto& msh = com_mod.msh[iM]; + int insd = nsd; + if (msh.lShl) insd = nsd-1; + if (msh.lFib) insd = 1; + + // Update pressure function space for Taylor-Hood type element + if (com_mod.msh[iM].nFs == 2) { + fs.nG = msh.fs[1].nG; + fs.eType = msh.fs[1].eType; + fs.lShpF = msh.fs[1].lShpF; + fs.eNoN = msh.fs[1].eNoN; + + fs.w.resize(fs.nG); + fs.N.resize(fs.eNoN,fs.nG); + fs.Nx.resize(nsd,fs.eNoN,fs.nG); + + if (fs.eType != ElementType::NRB) { + fs.w = msh.fs[1].w; + fs.N = msh.fs[1].N; + fs.Nx = msh.fs[1].Nx; + } + + } else { + fs.nG = msh.fs[0].nG; + fs.eType = msh.fs[0].eType; + fs.lShpF = msh.fs[0].lShpF; + fs.eNoN = msh.fs[0].eNoN; + + fs.w.resize(fs.nG); + fs.N.resize(fs.eNoN,fs.nG); + fs.Nx.resize(nsd,fs.eNoN,fs.nG); + + if (fs.eType != ElementType::NRB) { + fs.w = msh.fs[0].w; + fs.N = msh.fs[0].N; + fs.Nx = msh.fs[0].Nx; + } + } + int eNoN = fs.eNoN; + + Array xl(nsd,eNoN); + Array Nxi(insd,eNoN); + Array Nx(insd,eNoN); + Vector sl(eNoN); + Array tmps(nsd,insd); + Array tmp(nsd,nsd); + + for (int e = 0; e < msh.nEl; e++) { + // Updating the shape functions, if this is a NURB + // + // [TODO:DaveP] not implemented. + // + if (msh.eType == ElementType::NRB) { + //CALL NRBNNX(msh(iM), e) + fs.w = msh.w; + fs.N = msh.N; + fs.Nx = msh.Nx; + } + + int ibl = 0; + for (int a = 0; a < eNoN; a++) { + int Ac = msh.IEN(a,e); + xl.set_col(a, com_mod.x.col(Ac)); + + if (com_mod.mvMsh) { + for (int i = 0; i < nsd; i++) { + xl(i,a) += com_mod.Do(i+nsd+1,Ac); + } + } + sl(a) = s(0,Ac); + ibl = ibl + com_mod.iblank(Ac); + } + + if (ibl == eNoN) { + continue; + } + + double Jac = 0.0; + + for (int g = 0; g < fs.nG; g++) { + Nxi = fs.Nx.slice(g); + + if (g == 0 || !fs.lShpF) { + if (msh.lShl) { + Vector nV(nsd); + nn::gnns(nsd, eNoN, Nxi, xl, nV, tmps, tmps); + Jac = sqrt(utils::norm(nV)); + } else { + nn::gnn(eNoN, nsd, insd, Nxi, xl, Nx, Jac, tmp); + } + } + + if (utils::is_zero(Jac)) { + throw std::runtime_error("Jac < 0 for element: " + std::to_string(e) + ")"); + } + double sHat = 0.0; + + for (int a = 0; a < eNoN; a++) { + int Ac = msh.IEN(a,e); + sHat = sHat + sl(a)*fs.N(a,g); + } + result += fs.w(g) * Jac * sHat; + } + } + + if (com_mod.cm.seq()) { + return result; + } + result = com_mod.cm.reduce(cm_mod, result); + return result; +} +/////////////////////////////////////////////////////////////////////////// + /// @brief This routine integrated a scalar field over a particular domain. /// /// Note that 'l' and 'u' should be 0-based and are used to index into 's'. diff --git a/Code/Source/solver/all_fun.h b/Code/Source/solver/all_fun.h index 681646d80..c2dea3ba6 100644 --- a/Code/Source/solver/all_fun.h +++ b/Code/Source/solver/all_fun.h @@ -56,6 +56,8 @@ namespace all_fun { Array global(const ComMod& com_mod, const CmMod& cm_mod, const mshType& lM, const Array& U); + double integ(const ComMod& com_mod, const CmMod& cm_mod, int iM, const Array& s); + double integ(const ComMod& com_mod, const CmMod& cm_mod, int dId, const Array& s, int l, int u, bool pFlag=false); diff --git a/Code/Source/solver/consts.h b/Code/Source/solver/consts.h index 4d3b35c4c..fe4cc6136 100644 --- a/Code/Source/solver/consts.h +++ b/Code/Source/solver/consts.h @@ -128,7 +128,8 @@ enum class BoundaryConditionType bType_free = 19, // shell free bType_symm = 20, // shell symmetric bType_undefNeu = 21, // undeforming Neu - bType_RCR = 22 // RCR-Neu + bType_RCR = 22, // RCR-Neu + bType_Ris0D = 23, // RIS 0D }; // Define constants using smaller name and integer value (needed for bitwise operations). @@ -190,6 +191,9 @@ constexpr auto iBC_undefNeu = static_cast(BoundaryConditionType::bType_unde constexpr auto BC_ustd = BoundaryConditionType::bType_ustd; constexpr auto iBC_ustd = static_cast(BoundaryConditionType::bType_ustd); +constexpr auto BC_Ris0D = BoundaryConditionType::bType_Ris0D; +constexpr auto iBC_Ris0D = static_cast(BoundaryConditionType::bType_Ris0D); + //----------------------- // ConstitutiveModelType //----------------------- diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index a117ed918..a503f2b2d 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -172,6 +172,21 @@ void distribute(Simulation* simulation) dmsg << "wgt: " << wgt; #endif + // Need the ris flag here + cm.bcast(cm_mod, &com_mod.risFlag); + // Communicating RIS info + if (com_mod.risFlag) { + dist_ris(com_mod, cm_mod, cm); + } + // Need the uris flag here + cm.bcast(cm_mod, &com_mod.urisFlag); + if (com_mod.urisFlag) { + dist_uris(com_mod, cm_mod, cm); + } + + + int risProc = -1; + int jM = 0; for (int iM = 0; iM < nMsh; iM++) { #ifdef debug_distribute dmsg << " " << " "; @@ -184,6 +199,34 @@ void distribute(Simulation* simulation) #ifdef debug_distribute dmsg << "iWgt: " << iWgt; #endif + + if (com_mod.risFlag) { + for (int iProj = 0; iProj < com_mod.ris.nbrRIS; iProj++) { + // Find the adjacent mesh + if (com_mod.ris.lst(0,0,iProj) == iM) { + jM = com_mod.ris.lst(1,0,iProj); + } else if (com_mod.ris.lst(1,0,iProj) == iM) { + jM = com_mod.ris.lst(0,0,iProj); + } else { + continue; + } + // Right now just assign all nodes to the same proc + risProc = -1; + for (int e = 0; e < com_mod.msh[jM].gnEl; e++) { + if (com_mod.msh[jM].eRIS(e)) { + risProc = com_mod.msh[jM].partRIS(e); + break; + } + } + if (risProc != -1) { + for (int e = 0; e < com_mod.msh[iM].gnEl; e++) { + if (com_mod.msh[iM].eRIS(e)) { + com_mod.msh[iM].partRIS(e) = risProc; + } + } + } + } + } part_msh(simulation, iM, com_mod.msh[iM], gmtl, num_proc, iWgt); } @@ -288,6 +331,7 @@ void distribute(Simulation* simulation) dmsg << "Sending data read by master to slaves " << " ..."; #endif + if (!com_mod.resetSim) { cm.bcast(cm_mod, &com_mod.nsymd); cm.bcast(cm_mod, com_mod.stopTrigName); @@ -321,7 +365,12 @@ void distribute(Simulation* simulation) cm.bcast(cm_mod, &com_mod.sstEq); cm.bcast(cm_mod, &simulation->cep_mod.cepEq); - + cm.bcast(cm_mod, &com_mod.risFlag); + cm.bcast(cm_mod, &com_mod.ris0DFlag); + cm.bcast(cm_mod, &com_mod.urisFlag); + cm.bcast(cm_mod, &com_mod.urisActFlag); + cm.bcast(cm_mod, &com_mod.urisRes); + cm.bcast(cm_mod, &com_mod.urisResClose); cm.bcast(cm_mod, &com_mod.usePrecomp); if (com_mod.rmsh.isReqd) { auto& rmsh = com_mod.rmsh; @@ -551,6 +600,22 @@ void distribute(Simulation* simulation) } cm.bcast(cm_mod, &cplBC.initRCR); + + if (com_mod.risFlag) { + for (int iProj = 0; iProj < com_mod.ris.nbrRIS; iProj++) { + for (int i = 0; i < 2; i++) { + for (int j = 0; j < com_mod.risMapList[iProj].map.ncols(); j++) { + int tmp = gmtl(com_mod.grisMapList[iProj].map(i,j)); + if (tmp != -1) { + com_mod.risMapList[iProj].map(i,j) = com_mod.msh[i].lN(tmp); + } else { + com_mod.risMapList[iProj].map(i,j) = -1; + } + com_mod.grisMapList[iProj].map(i,j) = tmp; + } + } + } + } } @@ -569,6 +634,7 @@ void dist_bc(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bcType& lBc bool is_slave = cm.slv(cm_mod); cm.bcast(cm_mod, &lBc.cplBCptr); cm.bcast(cm_mod, &lBc.bType); + cm.bcast(cm_mod, &lBc.clsFlgRis); int nsd = com_mod.nsd; #ifdef debug_dist_bc @@ -800,6 +866,7 @@ void dist_bc(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bcType& lBc } } + //--------- // dist_bf //--------- @@ -937,6 +1004,466 @@ void dist_bf(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bfType& lBf } } +//--------- +// dist_ris +//--------- +/// @brief This rountine distributes data structures for RIS +void dist_ris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) +{ + using namespace consts; + + #define n_debug_dist_ris + #ifdef debug_dist_ris + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + int task_id = cm.idcm(); + auto& RIS = com_mod.ris; + + #ifdef debug_dist_ris + dmsg << "risFlag: " << com_mod.risFlag; + #endif + + Vector dims(3); + int lst_size = 0; + int maplist_size = 0; + Vector nStks; + int nProj = 0; + + // First find the RIS nodes between mesh pairs on the current processor + if (cm.mas(cm_mod)) { + nProj = RIS.nbrRIS; + dims(0) = RIS.lst.nrows(); + dims(1) = RIS.lst.ncols(); + dims(2) = RIS.lst.nslices(); + lst_size = RIS.lst.size(); + nStks.resize(nProj); + for (int iProj = 0; iProj < nProj; iProj++) { + nStks(iProj) = com_mod.risMapList[iProj].map.ncols(); + maplist_size += 2*nStks(iProj); + } + } + + cm.bcast(cm_mod, &nProj); + cm.bcast(cm_mod, dims); + cm.bcast(cm_mod, &lst_size); + if (cm.slv(cm_mod)) {nStks.resize(nProj);} + cm.bcast(cm_mod, nStks); + cm.bcast(cm_mod, &maplist_size); + + if (cm.slv(cm_mod)) { + com_mod.risMapList.resize(nProj); + com_mod.grisMapList.resize(nProj); + for (int iProj = 0; iProj < nProj; iProj++) { + com_mod.risMapList[iProj].map.resize(2, nStks(iProj)); + com_mod.grisMapList[iProj].map.resize(2, nStks(iProj)); + } + RIS.nbrIter.resize(nProj); + RIS.Res.resize(nProj); + RIS.clsFlg.resize(nProj); + RIS.meanP.resize(nProj,2); + RIS.meanFl.resize(nProj); + RIS.status.resize(nProj); + RIS.lst.resize(dims(0),dims(1),dims(2)); + } + dims.clear(); + + Vector flat_lst; + Vector flat_maplist; + Vector flat_gmaplist; + + if (cm.mas(cm_mod)) { + flat_lst.resize(lst_size); + int n = 0; + for (int i = 0; i < RIS.lst.nrows(); i++) { + for (int j = 0; j < RIS.lst.ncols(); j++) { + for (int k = 0; k < RIS.lst.nslices(); k++) { + flat_lst(n) = RIS.lst(i,j,k); + n += 1; + } + } + } + n = 0; + flat_maplist.resize(maplist_size); + flat_gmaplist.resize(maplist_size); + for (int iProj = 0; iProj < nProj; iProj++) { + for (int i = 0; i < 2; i++) { + for (int j = 0; j < nStks(iProj); j++) { + flat_maplist(n) = com_mod.risMapList[iProj].map(i,j); + flat_gmaplist(n) = com_mod.grisMapList[iProj].map(i,j); + n += 1; + } + } + } + } else { + flat_lst.resize(lst_size); + flat_maplist.resize(maplist_size); + flat_gmaplist.resize(maplist_size); + } + + cm.bcast(cm_mod, flat_lst); + cm.bcast(cm_mod, flat_maplist); + cm.bcast(cm_mod, flat_gmaplist); + cm.bcast(cm_mod, RIS.clsFlg); + cm.bcast(cm_mod, RIS.nbrIter); + cm.bcast(cm_mod, RIS.Res); + cm.bcast(cm_mod, &RIS.nbrRIS); + cm.bcast(cm_mod, RIS.meanP); + cm.bcast(cm_mod, RIS.meanFl); + cm.bcast(cm_mod, RIS.status); + + // Reshape the 1D array back to 3D after broadcasting + if (cm.slv(cm_mod)) { + int n = 0; + for (int i = 0; i < RIS.lst.nrows(); i++) { + for (int j = 0; j < RIS.lst.ncols(); j++) { + for (int k = 0; k < RIS.lst.nslices(); k++) { + RIS.lst(i,j,k) = flat_lst(n); + n += 1; + } + } + } + n = 0; + for (int iProj = 0; iProj < nProj; iProj++) { + for (int i = 0; i < 2; i++) { + for (int j = 0; j < nStks(iProj); j++) { + com_mod.risMapList[iProj].map(i,j) = flat_maplist(n); + com_mod.grisMapList[iProj].map(i,j) = flat_gmaplist(n); + n += 1; + } + } + } + } + flat_lst.clear(); + flat_maplist.clear(); + flat_gmaplist.clear(); + + // // [HZ] broadcasting nestd vector causes a problem in parallel + // for (int iProj = 0; iProj < nProj; iProj++) { + // cm.bcast(cm_mod, com_mod.risMapList[iProj].map); + // cm.bcast(cm_mod, com_mod.grisMapList[iProj].map); + // } + + for (int iM = 0; iM < com_mod.nMsh; iM++) { + cm.bcast(cm_mod, &com_mod.msh[iM].gnEl); + if (cm.slv(cm_mod)) { + com_mod.msh[iM].eRIS.resize(com_mod.msh[iM].gnEl); + com_mod.msh[iM].partRIS.resize(com_mod.msh[iM].gnEl); + } + cm.bcast(cm_mod, com_mod.msh[iM].eRIS); + cm.bcast(cm_mod, com_mod.msh[iM].partRIS); + } + nStks.clear(); +} + + +//--------- +// dist_uris +//--------- +/// @brief This routine distributes data structures for URIS +void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm) { + using namespace consts; + + #define n_debug_dist_uris + #ifdef debug_dist_uris + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + auto& uris = com_mod.uris; + + // We will copy the URIS valves to all processors + cm.bcast(cm_mod, &com_mod.nUris); + if (cm.slv(cm_mod)) { + uris.resize(com_mod.nUris); + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + uris[iUris].nrm.resize(com_mod.nsd); + } + } + + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + cm.bcast(cm_mod, uris[iUris].name); + cm.bcast(cm_mod, &uris[iUris].tnNo); + cm.bcast(cm_mod, &uris[iUris].nFa); + cm.bcast(cm_mod, &uris[iUris].sdf_default); + cm.bcast(cm_mod, &uris[iUris].sdf_deps); + cm.bcast(cm_mod, &uris[iUris].sdf_deps_close); + cm.bcast(cm_mod, &uris[iUris].clsFlg); + cm.bcast(cm_mod, &uris[iUris].cnt); + cm.bcast(cm_mod, &uris[iUris].scF); + cm.bcast(cm_mod, uris[iUris].nrm); + } + + std::vector> lM_gN_flat(com_mod.nUris); + std::vector> lM_lN_flat(com_mod.nUris); + std::vector> lM_IEN_flat(com_mod.nUris); + Vector lM_gN_flat_size(com_mod.nUris); + lM_gN_flat_size = 0; + Vector lM_lN_flat_size(com_mod.nUris); + lM_lN_flat_size = 0; + Vector lM_IEN_flat_size(com_mod.nUris); + lM_IEN_flat_size = 0; + int n0, n1, n2; + + if (cm.mas(cm_mod)) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + for (int iM = 0; iM < uris[iUris].nFa; iM++) { + auto& lM = uris[iUris].msh[iM]; + lM_gN_flat_size(iUris) += lM.gnNo; + lM_lN_flat_size(iUris) += uris[iUris].tnNo; + lM_IEN_flat_size(iUris) += lM.eNoN*lM.gnEl; + } + // std::cout << "lM_gN_flat_size: " << lM_gN_flat_size(iUris) << " for uris: " << iUris << std::endl; + // std::cout << "lM_lN_flat_size: " << lM_lN_flat_size(iUris) << " for uris: " << iUris << std::endl; + // std::cout << "lM_IEN_flat_size: " << lM_IEN_flat_size(iUris) << " for uris: " << iUris << std::endl; + } + } + + cm.bcast(cm_mod, lM_gN_flat_size); + cm.bcast(cm_mod, lM_lN_flat_size); + cm.bcast(cm_mod, lM_IEN_flat_size); + + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + lM_gN_flat[iUris].resize(lM_gN_flat_size(iUris)); + lM_lN_flat[iUris].resize(lM_lN_flat_size(iUris)); + lM_IEN_flat[iUris].resize(lM_IEN_flat_size(iUris)); + } + + if (cm.mas(cm_mod)) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + n0 = 0; + n1 = 0; + n2 = 0; + for (int iM = 0; iM < uris[iUris].nFa; iM++) { + auto& lM = uris[iUris].msh[iM]; + + for (int i = 0; i < lM.gnNo; i++) { + lM_gN_flat[iUris](n0) = lM.gN(i); + n0 += 1; + } + + for (int i = 0; i < uris[iUris].tnNo; i++) { + lM_lN_flat[iUris](n1) = lM.lN(i); + n1 += 1; + } + + for (int i = 0; i < lM.eNoN; i++) { + for (int j = 0; j < lM.gnEl; j++) { + lM_IEN_flat[iUris](n2) = lM.IEN(i,j); + n2 += 1; + } + } + + } + } + } + + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + cm.bcast(cm_mod, lM_gN_flat[iUris]); + cm.bcast(cm_mod, lM_lN_flat[iUris]); + cm.bcast(cm_mod, lM_IEN_flat[iUris]); + } + + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (cm.slv(cm_mod)) { + uris[iUris].msh.resize(uris[iUris].nFa); + } + for (int iM = 0; iM < uris[iUris].nFa; iM++) { + dist_uris_msh(com_mod, cm_mod, cm, uris[iUris].msh[iM], iUris); + } + } + + if (cm.slv(cm_mod)) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + n0 = 0; + n1 = 0; + n2 = 0; + for (int iM = 0; iM < uris[iUris].nFa; iM++) { + auto& lM = uris[iUris].msh[iM]; + + for (int i = 0; i < lM.gnNo; i++) { + lM.gN(i) = lM_gN_flat[iUris](n0); + n0 += 1; + } + + for (int i = 0; i < uris[iUris].tnNo; i++) { + lM.lN(i) = lM_lN_flat[iUris](n1); + n1 += 1; + } + + for (int i = 0; i < lM.eNoN; i++) { + for (int j = 0; j < lM.gnEl; j++) { + lM.IEN(i,j) = lM_IEN_flat[iUris](n2); + n2 += 1; + } + } + } + } + } + + Vector OpenN(com_mod.nUris); + Vector CloseN(com_mod.nUris); + std::vector> DxOpenFlat(com_mod.nUris); + std::vector> DxCloseFlat(com_mod.nUris); + Vector DxOpenFlatSize(com_mod.nUris); + Vector DxCloseFlatSize(com_mod.nUris); + int n; + if (cm.mas(cm_mod)) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + OpenN(iUris) = uris[iUris].DxOpen.nrows(); + CloseN(iUris) = uris[iUris].DxClose.nrows(); + DxOpenFlatSize(iUris) = uris[iUris].DxOpen.size(); + DxCloseFlatSize(iUris) = uris[iUris].DxClose.size(); + DxOpenFlat[iUris].resize(DxOpenFlatSize(iUris)); + DxCloseFlat[iUris].resize(DxCloseFlatSize(iUris)); + n = 0; + for (int i = 0; i < uris[iUris].DxOpen.nrows(); i++) { + for (int j = 0; j < uris[iUris].DxOpen.ncols(); j++) { + for (int k = 0; k < uris[iUris].DxOpen.nslices(); k++) { + DxOpenFlat[iUris](n) = uris[iUris].DxOpen(i,j,k); + n += 1; + } + } + } + n = 0; + for (int i = 0; i < uris[iUris].DxClose.nrows(); i++) { + for (int j = 0; j < uris[iUris].DxClose.ncols(); j++) { + for (int k = 0; k < uris[iUris].DxClose.nslices(); k++) { + DxCloseFlat[iUris](n) = uris[iUris].DxClose(i,j,k); + n += 1; + } + } + } + } + } + + cm.bcast(cm_mod, OpenN); + cm.bcast(cm_mod, CloseN); + cm.bcast(cm_mod, DxOpenFlatSize); + cm.bcast(cm_mod, DxCloseFlatSize); + + if (cm.slv(cm_mod)) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + uris[iUris].DxOpen.resize(OpenN(iUris), com_mod.nsd, uris[iUris].tnNo); + uris[iUris].DxClose.resize(CloseN(iUris), com_mod.nsd, uris[iUris].tnNo); + uris[iUris].x.resize(com_mod.nsd, uris[iUris].tnNo); + uris[iUris].Yd.resize(com_mod.nsd, uris[iUris].tnNo); + DxOpenFlat[iUris].resize(DxOpenFlatSize(iUris)); + DxCloseFlat[iUris].resize(DxCloseFlatSize(iUris)); + } + } + + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + cm.bcast(cm_mod, uris[iUris].x); + cm.bcast(cm_mod, uris[iUris].Yd); + cm.bcast(cm_mod, DxOpenFlat[iUris]); + cm.bcast(cm_mod, DxCloseFlat[iUris]); + } + + if (cm.slv(cm_mod)) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + n = 0; + for (int i = 0; i < uris[iUris].DxOpen.nrows(); i++) { + for (int j = 0; j < uris[iUris].DxOpen.ncols(); j++) { + for (int k = 0; k < uris[iUris].DxOpen.nslices(); k++) { + uris[iUris].DxOpen(i,j,k) = DxOpenFlat[iUris](n); + n += 1; + } + } + } + n = 0; + for (int i = 0; i < uris[iUris].DxClose.nrows(); i++) { + for (int j = 0; j < uris[iUris].DxClose.ncols(); j++) { + for (int k = 0; k < uris[iUris].DxClose.nslices(); k++) { + uris[iUris].DxClose(i,j,k) = DxCloseFlat[iUris](n); + n += 1; + } + } + } + } + } + +} + + +/// @brief This rountine distributes data structures for URIS mesh +void dist_uris_msh(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, mshType& lM, const int iUris) +{ + using namespace consts; + + #define n_debug_dist_uris_msh + #ifdef debug_dist_uris_msh + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + cm.bcast(cm_mod, lM.name); + cm.bcast(cm_mod, &lM.lShpF); + cm.bcast(cm_mod, &lM.lShl); + cm.bcast(cm_mod, &lM.lFib); + cm.bcast_enum(cm_mod, &lM.eType); + cm.bcast(cm_mod, &lM.eNoN); + cm.bcast(cm_mod, &lM.gnEl); + cm.bcast(cm_mod, &lM.gnNo); + cm.bcast(cm_mod, &lM.nFa); + cm.bcast(cm_mod, &lM.nFs); + cm.bcast(cm_mod, &lM.nFn); + cm.bcast(cm_mod, &lM.nG); + cm.bcast(cm_mod, &lM.vtkType); + cm.bcast(cm_mod, &lM.scF); + cm.bcast(cm_mod, &lM.dx); + // cm.bcast(cm_mod, &lM.fib.nFn); + // cm.bcast(cm_mod, &lM.fib.locNd); + // cm.bcast(cm_mod, &lM.fib.locEl); + // cm.bcast(cm_mod, &lM.fib.locGP); + + if (cm.slv(cm_mod)) { + lM.nNo = lM.gnNo; + lM.nEl = lM.gnEl; + lM.gN.resize(lM.nNo); + lM.lN.resize(com_mod.uris[iUris].tnNo); + lM.IEN.resize(lM.eNoN, lM.nEl); + lM.fa.resize(lM.nFa); + nn::select_ele(com_mod, lM); + // if (lM.fib.locNd) {lM.fib.fN.resize(nsd*lM.fib.nFn, 1, lM.nNo);} + // if (lM.fib.locEl) {lM.fib.fN.resize(nsd*lM.fib.nFn, 1, lM.nEl);} + // if (lM.fib.locGP) {lM.fib.fN.resize(nsd*lM.fib.nFn, lM.nG, lM.nEl);} + } + + // // [HZ] Broadcasting here causes problem, similar reason to the RIS map list + // cm.bcast(cm_mod, lM.gN); + // cm.bcast(cm_mod, lM.lN); + // cm.bcast(cm_mod, lM.IEN); + + // if (lM.fib.locNd || lM.fib.locEl || lM.fib.locGP) { + // cm.bcast(com_mod, &lM.fib.fN) + // } + + if (lM.eType == ElementType::NRB) { + cm.bcast(cm_mod, &lM.nSl); + int insd = com_mod.nsd; + if (lM.lShl) {insd = com_mod.nsd - 1;} + if (cm.slv(cm_mod)) { + lM.nW.resize(lM.gnNo); + lM.INN.resize(insd, lM.gnEl); + lM.bs.resize(insd); + } + cm.bcast(cm_mod, lM.nW); + cm.bcast(cm_mod, lM.INN); + for (int i = 0; i < insd; i++) { + cm.bcast(cm_mod, &lM.bs[i].n); + cm.bcast(cm_mod, &lM.bs[i].nG); + cm.bcast(cm_mod, &lM.bs[i].nEl); + cm.bcast(cm_mod, &lM.bs[i].nSl); + cm.bcast(cm_mod, &lM.bs[i].p); + if (cm.slv(cm_mod)) {lM.bs[i].xi.resize(lM.bs[i].n);} + cm.bcast(cm_mod, lM.bs[i].xi); + lM.bs[i].nNo = lM.bs[i].n - lM.bs[i].p - 1; + } + } + +} + void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std::vector& tMs, const Vector& gmtl, CepMod& cep_mod, eqType& lEq) { @@ -1652,6 +2179,7 @@ void part_msh(Simulation* simulation, int iM, mshType& lM, Vector& gmtl, in if (fp) fclose(fp); #endif + if (lM.eType == consts::ElementType::NRB) { part = cm.id(); @@ -1687,13 +2215,13 @@ void part_msh(Simulation* simulation, int iM, mshType& lM, Vector& gmtl, in #endif lM.IEN.resize(eNoN, nEl); + // Send lM.gIEN array to all processor's lM.IEN[] array of siize nEl*eNoN. // #ifdef dbg_part_msh dmsg << "sCount: " << sCount; dmsg << "disp: " << disp; #endif - MPI_Scatterv(lM.gIEN.data(), sCount.data(), disp.data(), cm_mod::mpint, lM.IEN.data(), nEl*eNoN, cm_mod::mpint, cm_mod.master, cm.com()); @@ -1704,6 +2232,7 @@ void part_msh(Simulation* simulation, int iM, mshType& lM, Vector& gmtl, in dmsg << "lM.IEN.size(): " << lM.IEN.size(); #endif + // The output of this process is "part" array which part(i) says // which processor element "i" belongs to // Doing partitioning, using ParMetis @@ -1781,10 +2310,35 @@ void part_msh(Simulation* simulation, int iM, mshType& lM, Vector& gmtl, in dmsg << " " << " "; dmsg << "Making the lM%IEN array " << " ..."; #endif - + if (cm.mas(cm_mod)) { sCount = 0; + int eRisProc = -1; for (int e = 0; e < lM.gnEl; e++) { + if (com_mod.risFlag) { + // If we found that this element needs to be shared + // since elements on the other mesh had been assigned + // a processor, then we change the gPart id + if (lM.partRIS(e) != -1) { + gPart[e] = lM.partRIS(e); + } else if (lM.eRIS(e)) { + // If this element is on a ris projection, + // we record the processor id so that next mesh + // will know that this element needs to be shared. + // Ideally, we might want to handle the case where + // elements next to the same projection are sent + // to different processors. Yet, this doesn't often + // happen and for simplicity, we + // force the processor id to be the same for all + // elements next to the same projection. + if (eRisProc == -1) { + eRisProc = gPart[e]; + } else { + gPart[e] = eRisProc; + } + lM.partRIS(e) = eRisProc; + } + } sCount[gPart[e]] = sCount[gPart[e]] + 1; } @@ -1853,6 +2407,9 @@ void part_msh(Simulation* simulation, int iM, mshType& lM, Vector& gmtl, in cm.bcast(cm_mod, &flag); cm.bcast(cm_mod, &fnFlag); cm.bcast(cm_mod, lM.eDist); + if (com_mod.risFlag) { + cm.bcast(cm_mod, lM.partRIS); + } nEl = lM.eDist[cm.id()+1] - lM.eDist[cm.id()]; #ifdef dbg_part_msh @@ -2088,5 +2645,4 @@ void part_msh(Simulation* simulation, int iM, mshType& lM, Vector& gmtl, in lM.Ys = all_fun::local(com_mod, cm_mod, cm, tmpYs); tmpYs.clear(); } -} - +} \ No newline at end of file diff --git a/Code/Source/solver/distribute.h b/Code/Source/solver/distribute.h index 33cfabd2c..97a84e799 100644 --- a/Code/Source/solver/distribute.h +++ b/Code/Source/solver/distribute.h @@ -40,6 +40,12 @@ void dist_bc(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bcType& lBc void dist_bf(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bfType& lBf); +void dist_ris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm); + +void dist_uris(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm); + +void dist_uris_msh(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, mshType& lM, const int iUris); + void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std::vector& tMs, const Vector& gmtl, CepMod& cep_mod, eqType& lEq); diff --git a/Code/Source/solver/fluid.cpp b/Code/Source/solver/fluid.cpp index a098a6670..c6b138c1e 100644 --- a/Code/Source/solver/fluid.cpp +++ b/Code/Source/solver/fluid.cpp @@ -36,6 +36,7 @@ #include "lhsa.h" #include "nn.h" #include "utils.h" +#include "ris.h" #include #include @@ -561,6 +562,8 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag // local tangent matrix (for a single element) Array3 lK(dof*dof,eNoN,eNoN); + double DDir = 0.0; + // Loop over all elements of mesh // int num_c = lM.nEl / 10; @@ -672,13 +675,43 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag dmsg << "w: " << w; #endif + // Plot the coordinates of the quad point in the current configuration + if (com_mod.urisFlag) { + Vector distSrf(com_mod.nUris); + distSrf = 0.0; + for (int a = 0; a < eNoN; a++) { + int Ac = lM.IEN(a,e); + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + distSrf(iUris) += fs[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); + } + } + + DDir = 0.0; + double DDirTmp = 0.0; + double sdf_deps_temp = 0.0; + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + if (com_mod.uris[iUris].clsFlg) { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + } else { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + } + if (distSrf(iUris) <= sdf_deps_temp) { + DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ + (2*sdf_deps_temp*sdf_deps_temp); + if (DDirTmp > DDir) {DDir = DDirTmp;} + } + } + + if (!com_mod.urisActFlag) {DDir = 0.0;} + } + // Compute momentum residual and tangent matrix. // if (nsd == 3) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); fluid_3d_m(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, - Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability); + Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); @@ -728,7 +761,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag if (nsd == 3) { auto N0 = fs[0].N.rcol(g); auto N1 = fs[1].N.rcol(g); - fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability); + fluid_3d_c(com_mod, vmsStab, fs[0].eNoN, fs[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, K_inverse_darcy_permeability, DDir); } else if (nsd == 2) { auto N0 = fs[0].N.rcol(g); @@ -739,6 +772,11 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag } // g: loop eq.linear_algebra->assemble(com_mod, eNoN, ptr, lK, lR); + if (com_mod.risFlag) { + if (!std::all_of(com_mod.ris.clsFlg.begin(), com_mod.ris.clsFlg.end(), [](bool v) { return v; })) { + ris::doassem_ris(com_mod, eNoN, ptr, lK, lR); + } + } } // e: loop @@ -1430,7 +1468,7 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability) + const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir) { #define n_debug_fluid3d_c #ifdef debug_fluid3d_c @@ -1456,6 +1494,14 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e const double ctM = 1.0; const double ctC = 36.0; + double Res; + if (!com_mod.urisFlag) { + Res = 0.0; + } else { + Res = com_mod.urisRes; + if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} + } + double rho = dmn.prop[PhysicalProperyType::fluid_density]; double f[3]; f[0] = dmn.prop[PhysicalProperyType::f_x]; @@ -1641,7 +1687,11 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e double kT = 4.0 * pow(ctM/dt,2.0); // If we consider the NSB model, we need to add an extra term inside the computation for the stab parameter - kT = kT + pow(K_inverse_darcy_permeability*mu/rho, 2.0); + kT = kT + pow(K_inverse_darcy_permeability*mu/rho, 2.0); + + // In case of unfitted RIS, compute the delta function at the quad point, + // add the additional value to the stabilization param + kT = kT + pow(Res*DDir, 2.0); double kU = u[0]*u[0]*Kxi(0,0) + u[1]*u[0]*Kxi(1,0) + u[2]*u[0]*Kxi(2,0) + u[0]*u[1]*Kxi(0,1) + u[1]*u[1]*Kxi(1,1) + u[2]*u[1]*Kxi(2,1) @@ -1664,13 +1714,25 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e rS[1] = mu_x[0]*es[0][1] + mu_x[1]*es[1][1] + mu_x[2]*es[2][1] + mu*d2u2[1]; rS[2] = mu_x[0]*es[0][2] + mu_x[1]*es[1][2] + mu_x[2]*es[2][2] + mu*d2u2[2]; - up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability*u[0]); - up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability*u[1]); - up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2]); + // up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability*u[0]); + // up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability*u[1]); + // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2]); + + up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability*u[0] + + (Res*DDir)*u[0]); + up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability*u[1] + + (Res*DDir)*u[1]); + up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability*u[2] + + (Res*DDir)*u[2]); for (int a = 0; a < eNoNw; a++) { double uNx = u[0]*Nwx(0,a) + u[1]*Nwx(1,a) + u[2]*Nwx(2,a); - T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a); + // T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a); + + T1 = -rho*uNx + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) + - mu*K_inverse_darcy_permeability*Nw(a) + - (Res*DDir)*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[1]*mu_g*esNx[0][a]; @@ -1738,7 +1800,7 @@ void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, - const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability) + const Array& bfl, Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir) { #define n_debug_fluid_3d_m #ifdef debug_fluid_3d_m @@ -1762,6 +1824,14 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e auto& dmn = eq.dmn[cDmn]; const double dt = com_mod.dt; + double Res; + if (!com_mod.urisFlag) { + Res = 0.0; + } else { + Res = com_mod.urisRes; + if (com_mod.uris[0].clsFlg) {Res = com_mod.urisResClose;} + } + double ctM = 1.0; double ctC = 36.0; @@ -1972,6 +2042,10 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // If we consider the NSB model, we need to add an extra term inside the computation for the stab parameter kT = kT + pow(K_inverse_darcy_permeability*mu/rho, 2.0); + // In case of unfitted RIS, compute the delta function at the quad point, + // add the additional value to the stabilization param + kT = kT + pow(Res*DDir, 2.0); + double kU = u[0]*u[0]*Kxi(0,0) + u[1]*u[0]*Kxi(1,0) + u[2]*u[0]*Kxi(2,0) + u[0]*u[1]*Kxi(0,1) + u[1]*u[1]*Kxi(1,1) + u[2]*u[1]*Kxi(2,1) + u[0]*u[2]*Kxi(0,2) + u[1]*u[2]*Kxi(1,2) + u[2]*u[2]*Kxi(2,2); @@ -2000,9 +2074,16 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e rS[2] = mu_x[0]*es[0][2] + mu_x[1]*es[1][2] + mu_x[2]*es[2][2] + mu*d2u2[2]; double up[3] = {}; - up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability * u[0]); - up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability * u[1]); - up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2]); + // up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability * u[0]); + // up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability * u[1]); + // up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2]); + + up[0] = -tauM*(rho*rV[0] + px[0] - rS[0] + mu*K_inverse_darcy_permeability * u[0] + + (Res*DDir)*u[0]); + up[1] = -tauM*(rho*rV[1] + px[1] - rS[1] + mu*K_inverse_darcy_permeability * u[1] + + (Res*DDir)*u[1]); + up[2] = -tauM*(rho*rV[2] + px[2] - rS[2] + mu*K_inverse_darcy_permeability * u[2] + + (Res*DDir)*u[2]); double tauC, tauB, pa; double eps = std::numeric_limits::epsilon(); @@ -2079,7 +2160,12 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e uaNx[a] = uNx[a]; } - T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a); + // T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) - mu*K_inverse_darcy_permeability*Nw(a); + + T1 = -rho*uNx[a] + mu*(Nwxx(0,a) + Nwxx(1,a) + Nwxx(2,a)) + + mu_x[0]*Nwx(0,a) + mu_x[1]*Nwx(1,a) + mu_x[2]*Nwx(2,a) + - mu*K_inverse_darcy_permeability*Nw(a) + - (Res*DDir)*Nw(a); updu[0][0][a] = mu_x[0]*Nwx(0,a) + d2u2[0]*mu_g*esNx[0][a] + T1; updu[1][0][a] = mu_x[1]*Nwx(0,a) + d2u2[1]*mu_g*esNx[0][a]; @@ -2114,7 +2200,9 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // dRm_a1/du_b1 double T2 = (mu + tauC)*rM[0][0] + esNx[0][a]*mu_g*esNx[0][b] - rho*tauM*uaNx[a]*updu[0][0][b]; lK(0,a,b) = lK(0,a,b) + wl*(T2 + T1); - lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); + // lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); + lK(0,a,b) = lK(0,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) + + (Res*DDir)*wl*Nw(b)*Nw(a); // dRm_a1/du_b2 T2 = mu*rM[1][0] + tauC*rM[0][1] + esNx[0][a]*mu_g*esNx[1][b] - rho*tauM*uaNx[a]*updu[1][0][b]; @@ -2131,7 +2219,9 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // dRm_a2/du_b2 T2 = (mu + tauC)*rM[1][1] + esNx[1][a]*mu_g*esNx[1][b] - rho*tauM*uaNx[a]*updu[1][1][b]; lK(5,a,b) = lK(5,a,b) + wl*(T2 + T1); - lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); + // lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); + lK(5,a,b) = lK(5,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) + + (Res*DDir)*wl*Nw(b)*Nw(a); // dRm_a2/du_b3 T2 = mu*rM[2][1] + tauC*rM[1][2] + esNx[1][a]*mu_g*esNx[2][b] - rho*tauM*uaNx[a]*updu[2][1][b]; @@ -2148,7 +2238,9 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // dRm_a3/du_b3; T2 = (mu + tauC)*rM[2][2] + esNx[2][a]*mu_g*esNx[2][b] - rho*tauM*uaNx[a]*updu[2][2][b]; lK(10,a,b) = lK(10,a,b) + wl*(T2 + T1); - lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); + // lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a); + lK(10,a,b) = lK(10,a,b) + mu*K_inverse_darcy_permeability*wl*Nw(b)*Nw(a) + + (Res*DDir)*wl*Nw(b)*Nw(a); //dmsg << "lK(10,a,b): " << lK(10,a,b); } } @@ -2173,10 +2265,14 @@ void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e // Residual contribution Birkman term // Local residue for (int a = 0; a < eNoNw; a++) { - lR(0,a) = lR(0,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[0]+up[0]); - lR(1,a) = lR(1,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[1]+up[1]); - lR(2,a) = lR(2,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[2]+up[2]); + lR(0,a) = lR(0,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[0]+up[0]) + + Res*DDir*w*Nw(a)*(u[0]+up[0]); + lR(1,a) = lR(1,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[1]+up[1]) + + Res*DDir*w*Nw(a)*(u[1]+up[1]); + lR(2,a) = lR(2,a) + mu*K_inverse_darcy_permeability*w*Nw(a)*(u[2]+up[2]) + + Res*DDir*w*Nw(a)*(u[2]+up[2]); } + } diff --git a/Code/Source/solver/fluid.h b/Code/Source/solver/fluid.h index 15c2c1e9b..be90a99ba 100644 --- a/Code/Source/solver/fluid.h +++ b/Code/Source/solver/fluid.h @@ -53,7 +53,7 @@ void construct_fluid(ComMod& com_mod, const mshType& lM, const Array& Ag void fluid_2d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability); + Array& lR, Array3& lK, double K_inverse_darcy_permeabilityx); void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, @@ -63,12 +63,12 @@ void fluid_2d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int e void fluid_3d_c(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability); + Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir=0.0); void fluid_3d_m(ComMod& com_mod, const int vmsFlag, const int eNoNw, const int eNoNq, const double w, const Array& Kxi, const Vector& Nw, const Vector& Nq, const Array& Nwx, const Array& Nqx, const Array& Nwxx, const Array& al, const Array& yl, const Array& bfl, - Array& lR, Array3& lK, double K_inverse_darcy_permeability); + Array& lR, Array3& lK, double K_inverse_darcy_permeability, double DDir=0.0); void get_viscosity(const ComMod& com_mod, const dmnType& lDmn, double& gamma, double& mu, double& mu_s, double& mu_x); diff --git a/Code/Source/solver/fsi.cpp b/Code/Source/solver/fsi.cpp index ae5e99ea3..f0c5e509a 100644 --- a/Code/Source/solver/fsi.cpp +++ b/Code/Source/solver/fsi.cpp @@ -38,6 +38,7 @@ #include "nn.h" #include "sv_struct.h" #include "utils.h" +#include "ris.h" #include #include @@ -103,6 +104,7 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar // double struct_3d_time = 0.0; double fluid_3d_time = 0.0; + double DDir = 0.0; for (int e = 0; e < lM.nEl; e++) { // setting globals @@ -210,6 +212,45 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar double w = fs_1[0].w(g) * Jac; + // Plot the coordinates of the quad point in the current configuration + if (com_mod.urisFlag) { + Vector distSrf(com_mod.nUris); + distSrf = 0.0; + for (int a = 0; a < eNoN; a++) { + int Ac = lM.IEN(a,e); + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + distSrf(iUris) += fs_1[0].N(a,g) * std::fabs(com_mod.uris[iUris].sdf(Ac)); + } + } + + DDir = 0.0; + double sdf_deps_temp = 0; + double DDirTmp = 0.0; + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + // if (distSrf(iUris) <= com_mod.uris[iUris].sdf_deps) { + // DDirTmp = (1 + cos(pi*distSrf(iUris)/com_mod.uris[iUris].sdf_deps))/ + // (2*com_mod.uris[iUris].sdf_deps*com_mod.uris[iUris].sdf_deps); + // if (DDirTmp > DDir) {DDir = DDirTmp;} + // } + + if (com_mod.uris[iUris].clsFlg) { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps_close; + } else { + sdf_deps_temp = com_mod.uris[iUris].sdf_deps; + } + if (distSrf(iUris) <= sdf_deps_temp) { + DDirTmp = (1 + cos(pi*distSrf(iUris)/sdf_deps_temp))/ + (2*sdf_deps_temp*sdf_deps_temp); + if (DDirTmp > DDir) {DDir = DDirTmp;} + } + } + + if (!com_mod.urisActFlag) {DDir = 0.0;} + + // std::cout << "===== DDir: " << DDir << std::endl; + } + + if (nsd == 3) { switch (cPhys) { case Equation_fluid: { @@ -217,7 +258,9 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar auto N1 = fs_1[1].N.col(g); // using zero permeability to use Navier-Stokes here, not Navier-Stokes-Brinkman - fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); + // fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); + fluid::fluid_3d_m(com_mod, vmsStab, fs_1[0].eNoN, fs_1[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir); + } break; case Equation_struct: { @@ -293,7 +336,8 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar auto N1 = fs_2[1].N.col(g); // using zero permeability to use Navier-Stokes here, not Navier-Stokes-Brinkman - fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); + //fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0); + fluid::fluid_3d_c(com_mod, vmsStab, fs_2[0].eNoN, fs_2[1].eNoN, w, ksix, N0, N1, Nwx, Nqx, Nwxx, al, yl, bfl, lR, lK, 0.0, DDir); } break; case Equation_ustruct: @@ -322,6 +366,11 @@ void construct_fsi(ComMod& com_mod, CepMod& cep_mod, const mshType& lM, const Ar eq.linear_algebra->assemble(com_mod, eNoN, ptr, lK, lR); + if (com_mod.risFlag) { + if (!std::all_of(com_mod.ris.clsFlg.begin(), com_mod.ris.clsFlg.end(), [](bool v) { return v; })) { + ris::doassem_ris(com_mod, eNoN, ptr, lK, lR); + } + } } // e: loop #ifdef debug_construct_fsi diff --git a/Code/Source/solver/initialize.cpp b/Code/Source/solver/initialize.cpp index a47e5b5a2..e871c4bda 100644 --- a/Code/Source/solver/initialize.cpp +++ b/Code/Source/solver/initialize.cpp @@ -95,6 +95,8 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< bool sstEq = com_mod.sstEq; bool pstEq = com_mod.pstEq; bool cepEq = cep_mod.cepEq; + bool risFlag = com_mod.risFlag; + bool urisFlag = com_mod.urisFlag; com_mod.timeP[0] = timeP[0]; com_mod.timeP[1] = timeP[1]; @@ -146,6 +148,21 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< bin_file.read((char*)Ad.data(), Ad.msize()); bin_file.read((char*)Xion.data(), Xion.msize()); bin_file.read((char*)cem.Ya.data(), cem.Ya.msize()); + } else if (risFlag) { + bin_file.read((char*)Ad.data(), Ad.msize()); + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} + } else if (urisFlag) { + bin_file.read((char*)Ad.data(), Ad.msize()); + Vector urisCnt(com_mod.nUris); + bin_file.read((char*)urisCnt.data(), urisCnt.msize()); + std::vector urisClsFlagChar(com_mod.nUris); + bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.nUris; i++) { + com_mod.uris[i].cnt = urisCnt(i); + com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} } else { bin_file.read((char*)Ad.data(), Ad.msize()); } @@ -156,6 +173,19 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< } else if (cepEq) { bin_file.read((char*)Xion.data(), Xion.msize()); bin_file.read((char*)cem.Ya.data(), cem.Ya.msize()); + } else if (risFlag) { + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} + } else if (urisFlag) { + Vector urisCnt(com_mod.nUris); + bin_file.read((char*)urisCnt.data(), urisCnt.msize()); + std::vector urisClsFlagChar(com_mod.nUris); + bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.nUris; i++) { + com_mod.uris[i].cnt = urisCnt(i); + com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} } else { //READ(fid,REC=cm.tF()) tStamp, cTS, time, timeP(1), eq.iNorm, cplBC.xo, Yo, Ao, Do } @@ -166,6 +196,19 @@ void init_from_bin(Simulation* simulation, const std::string& fName, std::array< } else { if (cepEq) { bin_file.read((char*)Xion.data(), Xion.msize()); + } else if (risFlag) { + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + bin_file.read(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + com_mod.ris.clsFlg[i] = clsFlagChar[i] ? 1 : 0;} + } else if (urisFlag) { + Vector urisCnt(com_mod.nUris); + bin_file.read((char*)urisCnt.data(), urisCnt.msize()); + std::vector urisClsFlagChar(com_mod.nUris); + bin_file.read(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); + for (int i = 0; i < com_mod.nUris; i++) { + com_mod.uris[i].cnt = urisCnt(i); + com_mod.uris[i].clsFlg = urisClsFlagChar[i] ? 1 : 0;} } else { //READ(fid,REC=cm.tF()) tStamp, cTS, time, timeP(1), eq.iNorm, cplBC.xo, Yo, Ao } @@ -509,6 +552,12 @@ void initialize(Simulation* simulation, Vector& timeP) i = i + cep_mod.nXion; if (cep_mod.cem.cpld) i = i + 1; } + if (com_mod.risFlag) { + i = i + com_mod.ris.nbrRIS; + } + if (com_mod.urisFlag) { + i = i + com_mod.nUris * 2; + } i = sizeof(int)*(1+com_mod.stamp.size()) + sizeof(double)*(2 + com_mod.nEq + com_mod.cplBC.nX + i*com_mod.tnNo); @@ -527,7 +576,25 @@ void initialize(Simulation* simulation, Vector& timeP) // [TODO:DaveP] this is not implemented. // - if (com_mod.shlEq) { + // if (com_mod.shlEq) { + // for (int iM = 0; iM < com_mod.nMsh; iM++) { + // auto& msh = com_mod.msh[iM]; + // if (msh.lShl) { + // if (msh.eType == ElementType::NRB) { + // msh.sbc.resize(msh.eNoN, msh.nEl); + // //ALLOCATE(msh(iM)%eIEN(0,0)) + // //ALLOCATE(msh(iM)%sbc(msh(iM)%eNoN,msh(iM)%nEl)) + // msh.sbc = 0; + // } else if (msh.eType == ElementType::TRI3) { + // baf_ini_ns::set_shl_xien(simulation, msh); + // //CALL SETSHLXIEN(msh(iM)) + // } + // } + // } + // } + + + if (com_mod.shlEq or com_mod.urisActFlag) { for (int iM = 0; iM < com_mod.nMsh; iM++) { auto& msh = com_mod.msh[iM]; if (msh.lShl) { @@ -544,6 +611,27 @@ void initialize(Simulation* simulation, Vector& timeP) } } + if (com_mod.risFlag) { + if (cm.mas(cm_mod)) { + std::cout << "Finally the map is: " << std::endl; + for (int iProj = 0; iProj < com_mod.ris.nbrRIS; iProj++) { + std::cout << " -p " << cm.idcm() << " -pj " << iProj << " RIS node: " << std::endl; + Vector map0 = com_mod.grisMapList[iProj].map.row(0); + for (int mapi = 0; mapi < map0.size(); mapi++) { + std::cout << map0[mapi] << "\t"; + } + std::cout << " " << std::endl; + std::cout << " -p " << cm.idcm() << " -pj " << iProj << " RIS node: " << std::endl; + Vector map1 = com_mod.grisMapList[iProj].map.row(1); + for (int mapi = 0; mapi < map1.size(); mapi++) { + std::cout << map1[mapi] << "\t"; + } + std::cout << " " << std::endl; + } + } + } + + // Initialize tensor operations // mat_fun::ten_init(nsd); diff --git a/Code/Source/solver/lhsa.cpp b/Code/Source/solver/lhsa.cpp index 608a57f28..a265776e5 100644 --- a/Code/Source/solver/lhsa.cpp +++ b/Code/Source/solver/lhsa.cpp @@ -177,6 +177,7 @@ void lhsa(Simulation* simulation, int& nnz) Array uInd(mnnzeic, tnNo); uInd = -1; int nMsh = com_mod.nMsh; + bool risFlag = com_mod.risFlag; for (int iM = 0; iM < nMsh; iM++) { auto& msh = com_mod.msh[iM]; @@ -190,6 +191,33 @@ void lhsa(Simulation* simulation, int& nnz) int colN = msh.IEN(b,e); add_col(tnNo, rowN, colN, mnnzeic, uInd); } + + // Add extra connections for cooresponding nodes in case of RIS + int jMRIS = 0; + if (risFlag) { + auto& RIS = com_mod.ris; + // If rowN is in the list of ris nodes + // Loop through all projections and find the one associated with the current mesh + for (int iProj = 0; iProj < RIS.nbrRIS; iProj++) { + if (RIS.lst(0,0,iProj) == iM) { + jMRIS = 1; + } else if (RIS.lst(1,0,iProj) == iM) { + jMRIS = 0; + } else { + continue; // Skip to the next iteration + } + std::array mapIdx; + utils::find_loc(com_mod.grisMapList[iProj].map, rowN, mapIdx); + if (mapIdx[0] != -1) { + int rowNR = com_mod.grisMapList[iProj].map(jMRIS,mapIdx[1]); + if (rowNR == -1) continue; + for (int b = 0; b < msh.eNoN; b++) { + int colN = msh.IEN(b,e); + add_col(tnNo, rowNR, colN, mnnzeic, uInd); + } + } + } + } } } } diff --git a/Code/Source/solver/main.cpp b/Code/Source/solver/main.cpp index dac1eda4f..e1fdac1fb 100644 --- a/Code/Source/solver/main.cpp +++ b/Code/Source/solver/main.cpp @@ -53,6 +53,8 @@ #include "txt.h" #include "ustruct.h" #include "vtk_xml.h" +#include "ris.h" +#include "uris.h" #include #include @@ -359,6 +361,8 @@ void iterate_solution(Simulation* simulation) set_bc::set_bc_dir(com_mod, An, Yn, Dn); + if (com_mod.urisFlag) {uris::uris_calc_sdf(com_mod);} + iterate_precomputed_time(simulation); // Inner loop for Newton iteration @@ -488,6 +492,14 @@ void iterate_solution(Simulation* simulation) set_bc::set_bc_dir_w(com_mod, Yg, Dg); + if (com_mod.risFlag) { + ris::ris_resbc(com_mod, Yg, Dg); + } + + if (com_mod.ris0DFlag) { + ris::ris0d_bc(com_mod, cm_mod, Yg, Dg); + } + // Apply contact model and add its contribution to residual // if (com_mod.iCntct) { @@ -610,7 +622,6 @@ void iterate_solution(Simulation* simulation) #endif break; } - output::output_result(simulation, com_mod.timeP, 2, iEqOld); inner_count += 1; @@ -634,6 +645,31 @@ void iterate_solution(Simulation* simulation) } */ + if (com_mod.risFlag) { + ris::ris_meanq(com_mod, cm_mod); + ris::ris_status(com_mod, cm_mod); + if (cm.mas(cm_mod)) { + std::cout << "Iteration: " << com_mod.cTS << std::endl; + for (int iProj = 0; iProj < com_mod.ris.nbrRIS; iProj++) { + std::cout << "Status for RIS projection: " << iProj << std::endl; + std::cout << " RIS iteration: " << com_mod.ris.nbrIter(iProj) << std::endl; + std::cout << " Is the valve close? " << com_mod.ris.clsFlg[iProj] << std::endl; + std::cout << " The status is: " << com_mod.ris.status[iProj] << std::endl; + } + } + + if (!std::all_of(com_mod.ris.status.begin(), com_mod.ris.status.end(), [](bool s) { return s; })) { + if (std::any_of(com_mod.ris.nbrIter.begin(), com_mod.ris.nbrIter.end(), [](int iter) { return iter <= 1; })) { + if (cm.mas(cm_mod)) { + std::cout << "Valve status just changed. Do not update" << std::endl; + } + } else { + ris::ris_updater(com_mod, cm_mod); + } + // goto label_11; + } + } + // Saving the TXT files containing average and fluxes (or ECGs) #ifdef debug_iterate_solution dmsg << "Saving the TXT files containing average and fluxes ..." << std::endl; @@ -744,6 +780,41 @@ void iterate_solution(Simulation* simulation) //CALL IB_OUTCPUT() } + // [HZ] Part related to RIS0D + if (cEq == 0 && com_mod.ris0DFlag) { + ris::ris0d_status(com_mod, cm_mod); + } + + // [HZ] Part related to unfitted RIS + // If the valve is active, look at the pressure difference + if (com_mod.urisFlag) { + for (int iUris = 0; iUris < com_mod.nUris; iUris++) { + com_mod.uris[iUris].cnt++; + if (com_mod.uris[iUris].clsFlg) { + uris::uris_meanp(com_mod, cm_mod, iUris); + // if (com_mod.uris[iUris].cnt == 1) { + // // GOTO 11 // The GOTO Statement in the Fortran code + // } + } else { + uris::uris_meanv(com_mod, cm_mod, iUris); + } + if (cm.mas(cm_mod)) { + std::cout << " URIS surface: " << com_mod.uris[iUris].name << ", count: " << com_mod.uris[iUris].cnt << std::endl; + } + } + + if (com_mod.mvMsh) { + uris::uris_update_disp(com_mod, cm_mod); + } + + if (cm.mas(cm_mod)) { + if (l2 && l3) { + uris::uris_write_vtus(com_mod); + } + } + } + // end RIS/URIS stuff + // Exiting outer loop if l1 if (l1) { break; @@ -823,8 +894,7 @@ int main(int argc, char *argv[]) dmsg << "Read files " << " ... "; #endif read_files(simulation, file_name); - - + // Distribute data to processors. #ifdef debug_main dmsg << "Distribute data to processors " << " ... "; @@ -881,5 +951,4 @@ int main(int argc, char *argv[]) } MPI_Finalize(); -} - +} \ No newline at end of file diff --git a/Code/Source/solver/output.cpp b/Code/Source/solver/output.cpp index 00b5184e2..1aeea765e 100644 --- a/Code/Source/solver/output.cpp +++ b/Code/Source/solver/output.cpp @@ -223,6 +223,8 @@ void write_restart(Simulation* simulation, std::array& timeP) const bool sstEq = com_mod.sstEq; const bool pstEq = com_mod.pstEq; const bool cepEq = cep_mod.cepEq; + const bool risFlag = com_mod.risFlag; + const bool urisFlag = com_mod.urisFlag; const auto& stFileName = com_mod.stFileName; auto& cplBC = com_mod.cplBC; @@ -309,6 +311,21 @@ void write_restart(Simulation* simulation, std::array& timeP) restart_file.write((char*)Ad.data(), Ad.msize()); restart_file.write((char*)Xion.data(), Xion.msize()); restart_file.write((char*)cem.Ya.data(), cem.Ya.msize()); + } else if (risFlag) { + restart_file.write((char*)Ad.data(), Ad.msize()); + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + clsFlagChar[i] = com_mod.ris.clsFlg[i] ? 1 : 0;} + restart_file.write(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + } else if (urisFlag) { + restart_file.write((char*)Ad.data(), Ad.msize()); + Vector urisCnt(com_mod.nUris); + std::vector urisClsFlagChar(com_mod.nUris); + for (int i = 0; i < com_mod.nUris; i++) { + urisCnt(i) = com_mod.uris[i].cnt; + urisClsFlagChar[i] = com_mod.uris[i].clsFlg ? 1 : 0;} + restart_file.write((char*)urisCnt.data(), urisCnt.msize()); + restart_file.write(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); } else { restart_file.write((char*)Ad.data(), Ad.msize()); } @@ -319,12 +336,40 @@ void write_restart(Simulation* simulation, std::array& timeP) } else if (cepEq) { restart_file.write((char*)Xion.data(), Xion.msize()); restart_file.write((char*)cem.Ya.data(), cem.Ya.msize()); + } else if (risFlag) { + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + clsFlagChar[i] = com_mod.ris.clsFlg[i] ? 1 : 0;} + restart_file.write(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + } else if (urisFlag) { + Vector urisCnt(com_mod.nUris); + std::vector urisClsFlagChar(com_mod.nUris); + for (int i = 0; i < com_mod.nUris; i++) { + urisCnt(i) = com_mod.uris[i].cnt; + urisClsFlagChar[i] = com_mod.uris[i].clsFlg ? 1 : 0;} + restart_file.write((char*)urisCnt.data(), urisCnt.msize()); + restart_file.write(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); + } else { + restart_file.write((char*)Dn.data(), Dn.msize()); } } } else { if (cepEq) { restart_file.write((char*)Xion.data(), Xion.msize()); + } else if (risFlag) { + std::vector clsFlagChar(com_mod.ris.clsFlg.size()); + for (int i = 0; i < com_mod.ris.clsFlg.size(); i++) { + clsFlagChar[i] = com_mod.ris.clsFlg[i] ? 1 : 0;} + restart_file.write(clsFlagChar.data(), clsFlagChar.size()*sizeof(char)); + } else if (urisFlag) { + Vector urisCnt(com_mod.nUris); + std::vector urisClsFlagChar(com_mod.nUris); + for (int i = 0; i < com_mod.nUris; i++) { + urisCnt(i) = com_mod.uris[i].cnt; + urisClsFlagChar[i] = com_mod.uris[i].clsFlg ? 1 : 0;} + restart_file.write((char*)urisCnt.data(), urisCnt.msize()); + restart_file.write(urisClsFlagChar.data(), urisClsFlagChar.size()*sizeof(char)); } else { //WRITE(fid, REC=myID) stamp, cTS, time, CPUT()-timeP(1), eq%iNorm, cplBC%xn, Yn, An } diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index b02f76078..3cd3f32a7 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -205,6 +205,13 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, } else if (std::set{"Coupled Momentum","CMM"}.count(bc_type)) { lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_CMM)); + } else if (std::set{"RIS0D", "RIS0D"}.count(bc_type)) { + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Ris0D)); + auto& com_mod = simulation->com_mod; + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Neu)); + com_mod.ris0DFlag = true; + lBc.resistance = bc_params->resistance.value(); + } else { throw std::runtime_error("[read_bc] Unknown boundary condition type '" + bc_type + "'."); } diff --git a/Code/Source/solver/read_msh.cpp b/Code/Source/solver/read_msh.cpp index 822b4aefa..54bb981b6 100644 --- a/Code/Source/solver/read_msh.cpp +++ b/Code/Source/solver/read_msh.cpp @@ -43,6 +43,8 @@ #include "vtk_xml.h" #include "vtk_xml_parser.h" +#include "uris.h" + #include #include #include @@ -1256,10 +1258,17 @@ void read_msh(Simulation* simulation) } } + com_mod.gtnNo = 0; + + // Examining the existance of a RIS surface and setting risMap. + // Reseting gtnNo and recounting nodes that are not duplicated + set_ris_projector(simulation); + + set_uris_meshes(simulation); + // Examining the existance of projection faces and setting %gN. // Reseting gtnNo and recounting nodes that are not duplicated // - com_mod.gtnNo = 0; utils::stackType avNds; set_projector(simulation, avNds); @@ -1686,6 +1695,362 @@ void read_msh(Simulation* simulation) } } + +void set_ris_projector(Simulation* simulation) +{ + #define n_debug_set_ris_projector + #ifdef debug_set_ris_projector + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); + dmsg.banner(); + #endif + + auto& com_mod = simulation->get_com_mod(); + auto& RIS = com_mod.ris; + int nPrj = simulation->parameters.RIS_projection_parameters.size(); + + if (nPrj > 0) { + com_mod.risFlag = true; + RIS.nbrRIS = nPrj; + RIS.lst.resize(2,2,nPrj); + RIS.lst = 0; + RIS.nbrIter.resize(nPrj); + RIS.nbrIter = 0; + RIS.Res.resize(nPrj); + RIS.Res = 0.0; + RIS.clsFlg.resize(nPrj, false); + RIS.meanP.resize(nPrj, 2); + RIS.meanP = 0.0; + RIS.meanFl.resize(nPrj); + RIS.meanFl = 0.0; + RIS.status.resize(nPrj, false); + + for (int iM = 0; iM < com_mod.nMsh; iM++) { + com_mod.msh[iM].eRIS.resize(com_mod.msh[iM].gnEl); + com_mod.msh[iM].eRIS = 0; // Allocate and initialize to false (0) + com_mod.msh[iM].partRIS.resize(com_mod.msh[iM].gnEl); + com_mod.msh[iM].partRIS = -1; // Allocate and initialize to -1 + } + } + + if (!com_mod.risFlag) return; + com_mod.risMapList.resize(nPrj); + com_mod.grisMapList.resize(nPrj); + + // Calculate an upper limit for the number required stacks + // + int nStk = 0; + for (auto& params : simulation->parameters.RIS_projection_parameters) { + auto ctmpi = params->name(); + int iM, iFa; + all_fun::find_face(com_mod.msh, ctmpi, iM, iFa); + nStk = nStk + com_mod.msh[iM].fa[iFa].nNo; + } + std::vector stk(nStk); + utils::stackType lPrj; + #ifdef debug_set_ris_projector + dmsg << "nStk: " << nStk; + #endif + + // Match the nodal coordinates for each projection face. + int iProj = 0; + for (auto& params : simulation->parameters.RIS_projection_parameters) { + int iM, iFa; + auto ctmpi = params->name(); + all_fun::find_face(com_mod.msh, ctmpi, iM, iFa); + // auto& face1 = com_mod.msh[iM].fa[iFa]; + #ifdef debug_set_ris_projector + dmsg << "iM: " << iM; + dmsg << "iFa: " << iFa; + dmsg << "face1.name: " << face1.name; + #endif + com_mod.risMapList[iProj].map.resize(2,nStk); + com_mod.grisMapList[iProj].map.resize(2,nStk); + com_mod.risMapList[iProj].map = 0; + com_mod.grisMapList[iProj].map = 0; + + int jM, jFa; + auto ctmpj = params->project_from_face(); + all_fun::find_face(com_mod.msh, ctmpj, jM, jFa); + // auto& face2 = com_mod.msh[jM].fa[jFa]; + #ifdef debug_set_ris_projector + dmsg << "jM: " << jM; + dmsg << "jFa: " << jFa; + dmsg << "face2.name: " << face2.name; + #endif + + double tol = params->projection_tolerance(); + com_mod.msh[jM].tol = tol; + double res = params->resistance(); + com_mod.msh[jM].res = res; + RIS.Res(iProj) = res; + + match_nodes(com_mod, com_mod.msh[iM].fa[iFa], com_mod.msh[jM].fa[jFa], + com_mod.msh[jM].tol, nStk, com_mod.risMapList[iProj].map); + + RIS.lst(0,0,iProj) = iM; + RIS.lst(1,0,iProj) = jM; + RIS.lst(0,1,iProj) = iFa; + RIS.lst(1,1,iProj) = jFa; + + // std::cout << "** Need to match face: " << jFa << " of " << std::endl; + // std::cout << " mesh " << jM << " into face " << iFa << " of mesh " << iM << std::endl; + // std::cout << "** The nbr of nodes to proj is " << nStk << std::endl; + // std::cout << " ** The res is " << com_mod.msh[jM].res << std::endl; + + for (int a = 0; a < com_mod.msh[iM].gnNo; a++) { + if (com_mod.msh[iM].gN[a] == -1) { + com_mod.msh[iM].gN[a] = com_mod.gtnNo; + com_mod.gtnNo = com_mod.gtnNo + 1; + } + } + for (int a = 0; a < com_mod.msh[jM].gnNo; a++) { + if (com_mod.msh[jM].gN[a] == -1) { + com_mod.msh[jM].gN[a] = com_mod.gtnNo; + com_mod.gtnNo = com_mod.gtnNo + 1; + } + } + iProj++; + } + + // Building the ris map between corresponding node with total enumeration + for (int iProj = 0; iProj < nPrj; iProj++) { + for (int i = 0; i < 2; i++) { + // std::cout << "iProj: " << iProj << " mesh: " << i << std::endl; + for (int j = 0; j < nStk; j++) { + if (com_mod.risMapList[iProj].map(i,j) != -1) { + com_mod.grisMapList[iProj].map(i,j) = + com_mod.msh[RIS.lst(i,0,iProj)].gN[com_mod.risMapList[iProj].map(i,j)]; + // std::cout << "local node " << com_mod.risMapList[iProj].map(i,j) + // << " global " << com_mod.grisMapList[iProj].map(i,j) << std::endl; + } + } + } + } + + // Mark elements of mesh that are connected to the RIS surface + for (int iProj = 0; iProj < RIS.nbrRIS; iProj++) { + for (int j = 0; j < 2; j++) { + int iM = RIS.lst(j,0,iProj); + for (int e = 0; e < com_mod.msh[iM].gnEl; e++) { + for (int a = 0; a < com_mod.msh[iM].eNoN; a++) { + int Ac = com_mod.msh[iM].gIEN(a,e); + Ac = com_mod.msh[iM].gN[Ac]; + std::array mapIdx; + utils::find_loc(com_mod.grisMapList[iProj].map, Ac, mapIdx); + if (mapIdx[0] != -1) { + com_mod.msh[iM].eRIS(e) = 1; + } + } + } + } + } +} + +void match_nodes(const ComMod& com_mod, const faceType& lFa, const faceType& pFa, + const double ptol, const int nNds, Array& map) +{ + #define n_debug_match_nodes + #ifdef debug_match_nodes + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "lFa.name: " << lFa.name; + dmsg << "pFa.name: " << pFa.name; + #endif + + int iM = lFa.iM; + int iSh = 0; + for (int i = 0; i < iM; i++) { + iSh = iSh + com_mod.msh[i].gnNo; + } + + int jM = pFa.iM; + int jSh = 0; + for (int j = 0; j < jM; j++) { + jSh = jSh + com_mod.msh[j].gnNo; + } + #ifdef debug_match_nodes + dmsg << "iM: " << iM; + dmsg << "iSh: " << iSh; + dmsg << "jM: " << jM; + dmsg << "jSh: " << jSh; + #endif + + // std::cout << " mesh iM " << iM << " nbr nodes " << com_mod.msh[iM].gnNo << std::endl; + // std::cout << " mesh iM " << iM << " nodes start at tot nbr " << iSh << std::endl; + // std::cout << " mesh jM " << jM << " nbr nodes " << com_mod.msh[jM].gnNo << std::endl; + // std::cout << " mesh jM " << jM << " nodes start at tot nbr " << jSh << std::endl; + + double tol; + double eps = std::numeric_limits::epsilon(); + + if (utils::is_zero(ptol)) { + tol = 1.e3 * eps; + } else { + tol = ptol; + } + + // We want to have approximately 1000 nodes in each block. So we + // calculate nBkd, which is the number of separate blockes in each + // direction, based on that. + // + int a = pFa.nNo; + int nBkd = static_cast( pow(a/1000.0, 0.333) + 0.5) ; + if (nBkd == 0) { + nBkd = 1; + } + int nsd = com_mod.nsd; + int nBk = pow(nBkd, nsd); + #ifdef debug_match_faces + dmsg << "a: " << a; + dmsg << "nBkd: " << nBkd; + dmsg << "nBk: " << nBk; + #endif + + Vector nodeBlk(a); + std::vector blk(nBk); + + // Find the extents of the domain and size of each block. + // + auto lfa_nodes = lFa.gN + iSh; + auto pfa_nodes = pFa.gN + jSh; + Vector xMin(com_mod.nsd), xMax(com_mod.nsd); + + for (int i = 0; i < nsd; i++) { + auto lfa_coords = com_mod.x.rows(i, lfa_nodes); + auto pfa_coords = com_mod.x.rows(i, pfa_nodes); + xMin[i] = std::min(lfa_coords.min(), pfa_coords.min()); + xMax[i] = std::max(lfa_coords.max(), pfa_coords.max()); + + if (xMin[i] < 0.0) { + xMin[i] = xMin[i]*(1.0+eps); + } else { + xMin[i] = xMin[i]*(1.0-eps); + } + + if (xMax[i] < 0.0) { + xMax[i] = xMax[i]*(1.0-eps); + } else { + xMax[i] = xMax[i]*(1.0+eps); + } + } + + auto dx = (xMax - xMin) / static_cast(nBkd); + std::vector nFlt(nsd); + + for (int i = 0; i < nsd; i++) { + if (utils::is_zero(dx[i])) { + nFlt[i] = false; + } else { + nFlt[i] = true; + } + } + + // std::cout << " xMax " << xMax << std::endl; + // std::cout << " xMin " << xMin << std::endl; + // std::cout << " dx " << dx << std::endl; + + // Find an estimation for size of each block + // + for (int a = 0; a < pFa.nNo; a++) { + int Ac = pFa.gN[a] + jSh; + auto coord = com_mod.x.col(Ac); + int iBk = find_blk(nsd, nBkd, nFlt, xMin, dx, coord); + nodeBlk[a] = iBk; + blk[iBk].n = blk[iBk].n + 1; + } + + for (int iBk = 0; iBk < nBk; iBk++) { + blk[iBk].gN = Vector(blk[iBk].n); + blk[iBk].n = 0; + } + + for (int a = 0; a < pFa.nNo; a++) { + int Ac = pFa.gN[a]; + int iBk = nodeBlk[a]; + blk[iBk].gN(blk[iBk].n) = Ac; + blk[iBk].n = blk[iBk].n + 1; + } + + // Doing the calculation for every single node on this face. + // + int cnt = 0; + + for (int a = 0; a < lFa.nNo; a++) { + int Ac = lFa.gN[a]; + auto coord = com_mod.x.col(Ac+iSh); + int iBk = find_blk(nsd, nBkd, nFlt, xMin, dx, coord); + + // Check all nodes on the other face. + auto minS = std::numeric_limits::max(); + int i; + for (int b = 0; b < blk[iBk].n; b++) { + int Bc = blk[iBk].gN[b]; + if ((iM == jM) && (Ac == Bc)) { + continue; + } + + auto diff = com_mod.x.col(Bc+jSh) - com_mod.x.col(Ac+iSh); + double ds = sqrt(diff*diff); + + if (ds < minS) { + minS = ds; + i = Bc; + } + } + + int Bc = i; + + if (tol < 0.0) { + // std::cout << "adding connection (/Ac,Bc/)) = (" << Ac << ", " << Bc << ")" << std::endl; + map(0,cnt) = Ac; + map(1,cnt) = Bc; + cnt = cnt + 1; + } else if (minS < tol) { + // std::cout << "adding connection (/Ac,Bc/)) = (" << Ac << ", " << Bc << ")" << std::endl; + map(0,cnt) = Ac; + map(1,cnt) = Bc; + cnt = cnt + 1; + } + } + + #ifdef debug_match_faces + dmsg << "cnt: " << cnt; + dmsg << "lFa.nNo: " << lFa.nNo; + #endif + + if (cnt != lFa.nNo) { + throw std::runtime_error("Failed to find matching nodes between faces " + lFa.name + " and " + pFa.name + "."); + } + + std::cout << " Finally the map0 is: " << std::endl; + Vector map0 = map.row(0); + for (int i = 0; i < map0.size(); i++) { + std::cout << map(0,i) << " " << "\t"; + } + std::cout << std::endl; + std::cout << " Finally the map1 is: " << std::endl; + Vector map1 = map.row(1); + for (int i = 0; i < map1.size(); i++) { + std::cout << map(1,i) << " " << "\t"; + } + std::cout << std::endl; +} + +void set_uris_meshes(Simulation* simulation) +{ + #define n_debug_set_uris_meshes + #ifdef debug_set_uris_meshes + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); + dmsg.banner(); + #endif + + int nUris = simulation->parameters.URIS_mesh_parameters.size(); + + if (nUris > 0) { + uris::uris_read_msh(simulation); + } +} + /// @brief Read domain from a dat file /// /// Reproduces the Fortran 'SETDMNIDFF' subroutine. @@ -1894,5 +2259,4 @@ void set_projector(Simulation* simulation, utils::stackType& avNds) } } -}; - +}; \ No newline at end of file diff --git a/Code/Source/solver/read_msh.h b/Code/Source/solver/read_msh.h index 0712da0d3..95ef6fddf 100644 --- a/Code/Source/solver/read_msh.h +++ b/Code/Source/solver/read_msh.h @@ -72,6 +72,8 @@ namespace read_msh_ns { void load_var_ini(Simulation* simulation, const ComMod& com_mod); void match_faces(const ComMod& com_mod, const faceType& face1, const faceType& face2, const double tol, utils::stackType& lPrj); + void match_nodes(const ComMod& com_mod, const faceType& lFa, const faceType& pFa, + const double ptol, const int nNds, Array& map); void read_fib_nff(Simulation* simulation, mshType& mesh, const std::string& fName, const std::string& kwrd, const int idx); void read_msh(Simulation* simulation); @@ -79,8 +81,9 @@ namespace read_msh_ns { void set_dmn_id_ff(Simulation* simulation, mshType& mesh, const std::string& file_name); void set_dmn_id_vtk(Simulation* simulation, mshType& mesh, const std::string& file_name, const std::string& kwrd); void set_projector(Simulation* simulation, utils::stackType& avNds); - - + void set_ris_projector(Simulation* simulation); + void set_uris_meshes(Simulation* simulation); + }; #endif diff --git a/Code/Source/solver/ris.cpp b/Code/Source/solver/ris.cpp new file mode 100644 index 000000000..65af9631c --- /dev/null +++ b/Code/Source/solver/ris.cpp @@ -0,0 +1,570 @@ +/* Copyright (c) Stanford University, The Regents of the University of California, and others. + * + * All Rights Reserved. + * + * See Copyright-SimVascular.txt for additional details. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject + * to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS + * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER + * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include "ris.h" + +#include "all_fun.h" +#include "lhsa.h" +#include "mat_fun.h" +#include "nn.h" +#include "utils.h" +#include "set_bc.h" + +namespace ris { + +/// @brief This subroutine computes the mean pressure and flux on the ris surface +void ris_meanq(ComMod& com_mod, CmMod& cm_mod) +{ + #define n_debug_ris_meanq + #ifdef debug_ris_meanq + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "risFlag: " << com_mod.risFlag; + #endif + + using namespace consts; + + auto& cm = com_mod.cm; + auto& eq = com_mod.eq; + auto& RIS = com_mod.ris; + auto& msh = com_mod.msh; + + const int nsd = com_mod.nsd; + const int cEq = com_mod.cEq; + + auto& An = com_mod.An; + auto& Ad = com_mod.Ad; + auto& Dn = com_mod.Dn; + auto& Yn = com_mod.Yn; + + Array tmpV(maxNSD, com_mod.tnNo); + + // Get the number of projections from RIS + int nPrj = RIS.nbrRIS; + + int iEq = 0; + int m = 1; + int s = eq[iEq].s + nsd; + int e = s + m - 1; + + for (int j = 0; j < Yn.ncols(); j++) { + tmpV(0,j) = Yn(s,j); + } + + for (int iProj = 0; iProj < nPrj; iProj++){ + RIS.meanFl(iProj) = 0.0; + for (int i = 0; i < 2; i++) { // Assuming two meshes + RIS.meanP(iProj, i) = 0.0; + int iM = RIS.lst(i,0,iProj); + int iFa = RIS.lst(i,1,iProj); + double tmp = msh[iM].fa[iFa].area; + RIS.meanP(iProj,i) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0)/tmp; + } + } + + // For the velocity + m = nsd; + s = eq[iEq].s; + e = s + m - 1; + + for (int iProj = 0; iProj < nPrj; iProj++) { + // tmpV[0:m,:] = Yn[s:e,:]; + for (int i = 0; i < m; i++) { + for (int j = 0; j < Yn.ncols(); j++) { + tmpV(i,j) = Yn(s+i,j); + } + } + int iM = RIS.lst(0,0,iProj); + int iFa = RIS.lst(0,1,iProj); + RIS.meanFl(iProj) = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, m-1); + + if (cm.mas(cm_mod)) { + std::cout << "For RIS projection: " << iProj << std::endl; + std::cout << " The average pressure is: " << RIS.meanP(iProj,0) << ", " + << RIS.meanP(iProj,1) << std::endl; + std::cout << " The average flow is: " << RIS.meanFl(iProj) << std::endl; + } + } +} + +/// @brief Weak treatment of RIS resistance boundary conditions +void ris_resbc(ComMod& com_mod, const Array& Yg, const Array& Dg) +{ + using namespace consts; + #define n_debug_ris_resbc + #ifdef debug_ris_resbc + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "risFlag: " << com_mod.risFlag; + #endif + + // auto& eq = com_mod.eq[cEq]; + auto& eq = com_mod.eq; + auto& RIS = com_mod.ris; + auto& msh = com_mod.msh; + + const int nsd = com_mod.nsd; + const int cEq = com_mod.cEq; + const int cDmn = com_mod.cDmn; + int nPrj = RIS.nbrRIS; + + bcType lBc; + + for (int iProj = 0; iProj < nPrj; iProj++) { + if (!RIS.clsFlg[iProj]) {continue;} + + // Weak Dirichlet BC for fluid/FSI equations + lBc.weakDir = true; + lBc.tauB = RIS.Res[iProj]; + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Dir)); + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_std)); + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_flat)); + + lBc.eDrn.resize(nsd); + lBc.eDrn = 0; + + for (int i = 0; i < 2; i++) { // Always two meshes for one projection + int iM = RIS.lst(i,0,iProj); + int iFa = RIS.lst(i,1,iProj); + lBc.gx.resize(msh[iM].fa[iFa].nNo); + lBc.gx = 0.; //1.0; + auto& cPhys = eq[cEq].dmn[cDmn].phys; + + double Yg_max, Yg_min, Yg_sum; + double Dg_max, Dg_min, Dg_sum; + std::array fs; + + if (cPhys == EquationType::phys_fluid) { + // Build the correct BC + set_bc::set_bc_dir_wl(com_mod, lBc, msh[iM], msh[iM].fa[iFa], Yg, Dg); + } + lBc.gx.clear(); + } + lBc.eDrn.clear(); + } +} + + +void setbc_ris(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, + const Array& Yg, const Array& Dg) +{ + // [HZ] looks not needed in the current implementation +} + + +/// @brief This subroutine updates the resistance and activation flag for the +/// closed and open configurations of the RIS surfaces +void ris_updater(ComMod& com_mod, CmMod& cm_mod) +{ + #define n_debug_ris_updater + #ifdef debug_ris_updater + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "risFlag: " << com_mod.risFlag; + #endif + + // auto& eq = com_mod.eq[cEq]; + auto& cm = com_mod.cm; + auto& eq = com_mod.eq; + auto& RIS = com_mod.ris; + + int nPrj = RIS.nbrRIS; + + for (int iProj = 0; iProj < nPrj; iProj++) { + // The valve is closed check if it should open + if (RIS.clsFlg[iProj]) { + // OPENING CONDITION: Check condition on the pressure difference + if (RIS.meanP(iProj,0) > RIS.meanP(iProj,1)) { + RIS.clsFlg[iProj] = false; + if (cm.mas(cm_mod)) { + std::cout << "RIS Proj " << iProj << ": Going from close to open." << std::endl; + } + RIS.nbrIter(iProj) = 0; + // I needed to update the state variables when the valve + // goes from close to open to prevent the valve goes back + // to close at the next iteration. This is needed only for + // close to open and cannot be used for open to close. + com_mod.Ao = com_mod.An; + com_mod.Yo = com_mod.Yn; + if (com_mod.dFlag) {com_mod.Do = com_mod.Dn;} + com_mod.cplBC.xo = com_mod.cplBC.xn; + } + } else { + // The valve is open, check if it should close. + // CLOSING CONDITION: Check existence of a backflow + if (RIS.meanFl(iProj) < 0.) { + RIS.clsFlg[iProj] = true; + if (cm.mas(cm_mod)) { + std::cout << "RIS Proj " << iProj << ": Going from open to close." << std::endl; + } + RIS.nbrIter(iProj) = 0; + } + } + } +} + +/// @brief This subroutine will check the valve status if it is admissible +/// or not, if not admissible we recompute the iteration until it will be +void ris_status(ComMod& com_mod, CmMod& cm_mod) +{ + using namespace consts; + + #define n_debug_ris_status + #ifdef debug_ris_status + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "risFlag: " << com_mod.risFlag; + #endif + + // auto& eq = com_mod.eq[cEq]; + auto& cm = com_mod.cm; + auto& eq = com_mod.eq; + auto& RIS = com_mod.ris; + + int nPrj = RIS.nbrRIS; + + for (int iProj = 0; iProj < nPrj; iProj++) { + RIS.nbrIter(iProj) += 1; + RIS.status[iProj] = true; + // If the valve is closed, check the pressure difference, + // if the pressure difference is negative the valve should be open + // the status is then not admissible + if (RIS.clsFlg[iProj]) { + if (RIS.meanP(iProj,0) > RIS.meanP(iProj,1)) { + if (cm.mas(cm_mod)) { + std::cout << "RIS Proj " << iProj << ": **** Not admissible, it should be open ****" << std::endl; + } + RIS.status[iProj] = false; + } + } else { + // If the valve is open, chech the flow, + // if the flow is negative the valve should be closed + // the status is then not admissible + if (RIS.meanFl(iProj) < 0.) { + if (cm.mas(cm_mod)) { + std::cout << "RIS Proj " << iProj << ": **** Not admissible, it should be closed ****" << std::endl; + } + RIS.status[iProj] = false; + } + } + } +} + +/// @brief This subroutine assembels the element stiffness matrix into the +/// global stiffness matrix (Val sparse matrix formatted as a vector) +void doassem_ris(ComMod& com_mod, const int d, const Vector& eqN, + const Array3& lK, const Array& lR) +{ + using namespace consts; + + #define n_debug_doassem_ris + #ifdef debug_doassem_ris + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "risFlag: " << com_mod.risFlag; + #endif + + auto& cm = com_mod.cm; + auto& eq = com_mod.eq; + auto& RIS = com_mod.ris; + + const int nsd = com_mod.nsd; + int nPrj = RIS.nbrRIS; + int rowNadj = 0; + double val_sum; + + for (int iProj = 0; iProj < nPrj; iProj++) { + if (RIS.clsFlg[iProj]) {continue;} + + for (int a = 0; a < d; a++) { + int rowN = eqN(a); + if (rowN == -1) {continue;} + + std::array mapIdx; + utils::find_loc(com_mod.grisMapList[iProj].map, rowN, mapIdx); + + if (mapIdx[0] == -1) {continue;} + + for (int jM = 0; jM < 2; jM++) { + if (jM == mapIdx[0]) {continue;} + if (com_mod.grisMapList[iProj].map(jM, mapIdx[1]) == -1) {continue;} + rowNadj = com_mod.grisMapList[iProj].map(jM, mapIdx[1]); + } + + if (rowNadj == -1) {continue;} + + for (int i = 0; i < com_mod.R.nrows(); i++) { + com_mod.R(i,rowNadj) = com_mod.R(i,rowNadj) + lR(i,a); + } + + for (int b = 0; b < d; b++) { + int colN = eqN(b); + std::array mapIdxC; + utils::find_loc(com_mod.grisMapList[iProj].map, colN, mapIdxC); + + if (mapIdxC[0] != -1) { + for (int jM = 0; jM < 2; jM++) { + if (jM == mapIdxC[0]) {continue;} + if (com_mod.grisMapList[iProj].map(jM, mapIdxC[1]) == -1) {continue;} + colN = com_mod.grisMapList[iProj].map(jM, mapIdxC[1]); + } + } //else { + // continue; + // } + + if (colN == -1) {continue;} + int left = com_mod.rowPtr(rowNadj); + int right = com_mod.rowPtr(rowNadj+1); + int ptr = (right + left) / 2; + + while (colN != com_mod.colPtr(ptr)) { + if (colN > com_mod.colPtr(ptr)) { + left = ptr; + } else { + right = ptr; + } + ptr = (right + left) / 2; + } + + for (int i = 0; i < com_mod.Val.nrows(); i++) { + com_mod.Val(i,ptr) = com_mod.Val(i,ptr) + lK(i,a,b); + } + } + } + } +} + +/// @brief This subroutine assembels the element stiffness matrix into the +/// global stiffness matrix (Val sparse matrix formatted as a vector) +void doassem_velris(ComMod& com_mod, const int d, const Array& eqN, + const Array3& lK, const Array& lR) +{ + // [HZ] looks not needed in the current implementation +} + +void clean_r_ris(ComMod& com_mod) +{ + // [HZ] looks not needed in the current implementation +} + +void setbcdir_ris(ComMod& com_mod, Array& lA, Array& lY, Array& lD) +{ + // [HZ] looks not needed in the current implementation +} + +/// RIS0D code +void ris0d_bc(ComMod& com_mod, CmMod& cm_mod, const Array& Yg, const Array& Dg) +{ + using namespace consts; + + #define n_debug_ris0d_bc + #ifdef debug_ris0d_bc + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "ris0DFlag: " << com_mod.ris0DFlag; + #endif + + auto& cm = com_mod.cm; + auto& eq = com_mod.eq; + auto& msh = com_mod.msh; + + const int nsd = com_mod.nsd; + const int cEq = com_mod.cEq; + + bcType lBc; + + for (int iBc = 0; iBc < eq[cEq].nBc; iBc++) { + int iFa = eq[cEq].bc[iBc].iFa; + int iM = eq[cEq].bc[iBc].iM; + + if (!utils::btest(eq[cEq].bc[iBc].bType, iBC_Ris0D)) {continue;} + + if (eq[cEq].bc[iBc].clsFlgRis == 1) { + // Weak Dirichlet BC for fluid/FSI equations + lBc.weakDir = true; + lBc.tauB = eq[cEq].bc[iBc].resistance; + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_Dir)); + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_std)); + lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_flat)); + + lBc.eDrn.resize(nsd); + lBc.eDrn = 0; + + // Apply bc Dir + lBc.gx.resize(msh[iM].fa[iFa].nNo); + lBc.gx = 1.0; + set_bc::set_bc_dir_wl(com_mod, lBc, msh[iM], msh[iM].fa[iFa], Yg, Dg); + lBc.gx.clear(); + lBc.eDrn.clear(); + } else { + // Apply Neu bc + set_bc::set_bc_neu_l(com_mod, cm_mod, eq[cEq].bc[iBc], msh[iM].fa[iFa], Yg, Dg); + } + + } + +} + +void ris0d_status(ComMod& com_mod, CmMod& cm_mod)//, const Array& Yg, const Array& Dg) +{ + using namespace consts; + + #define n_debug_status + #ifdef debug_status + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "ris0DFlag: " << com_mod.ris0DFlag; + #endif + + auto& cm = com_mod.cm; + auto& eq = com_mod.eq; + auto& msh = com_mod.msh; + + const int nsd = com_mod.nsd; + const int cEq = com_mod.cEq; + + auto& An = com_mod.An; + auto& Ad = com_mod.Ad; + auto& Dn = com_mod.Dn; + auto& Yn = com_mod.Yn; + + bcType lBc; + faceType lFa; + + double meanP = 0.0; + double meanFl = 0.0; + double tmp, tmp_new; + Vector sA; + Array tmpV; + + for (int iBc = 0; iBc < eq[cEq].nBc; iBc++) { + int iFa = eq[cEq].bc[iBc].iFa; + int iM = eq[cEq].bc[iBc].iM; + + if (!utils::btest(eq[cEq].bc[iBc].bType, iBC_Ris0D)) {continue;} + + tmpV.resize(maxNSD, com_mod.tnNo); + tmpV = 0.0; + + // Compute mean Q and pressure difference + int m = 1; + int s = eq[cEq].s + nsd; + // int e = s + m - 1; + + for (int j = 0; j < Yn.ncols(); j++) { + tmpV(0,j) = Yn(s,j); + } + + tmp = msh[iM].fa[iFa].area; + sA.resize(com_mod.tnNo); + sA = 1.0; + lFa = msh[iM].fa[iFa]; + // such update may be not correct + tmp_new = all_fun::integ(com_mod, cm_mod, lFa, sA); + meanP = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, m-1)/tmp_new; + + // For the velocity + m = nsd; + s = eq[cEq].s; + // e = s + m - 1; + tmpV.resize(maxNSD,com_mod.tnNo); + tmpV = 0.0; + + for (int i = 0; i < m; i++) { + for (int j = 0; j < Yn.ncols(); j++) { + tmpV(i,j) = Yn(s+i,j); + } + } + + meanFl = all_fun::integ(com_mod, cm_mod, msh[iM].fa[iFa], tmpV, 0, m-1); + + std::cout << "The average pressure is: " << meanP << std::endl; + std::cout << "The pressure from 0D is: " << eq[cEq].bc[iBc].g << std::endl; + std::cout << "The average flow is: " << meanFl << std::endl; + + com_mod.RisnbrIter = com_mod.RisnbrIter + 1; + + if (com_mod.RisnbrIter < 25 && com_mod.cTS > 0) { + if (!cm.seq()) {cm.bcast(cm_mod, &com_mod.RisnbrIter);} + return ; + } + + // Update RES + // Update the resistance - determine the configuration + // The valve is closed check if it should open + if (eq[cEq].bc[iBc].clsFlgRis == 1) { + // OPENING CONDITION: Check condition on the pressure difference + if (eq[cEq].bc[iBc].g < meanP) { + eq[cEq].bc[iBc].clsFlgRis = 0; + if (!cm.seq()) { + cm.bcast(cm_mod, &eq[cEq].bc[iBc].clsFlgRis); + } + std::cout << "!!! -- Going from close to open " << std::endl; + com_mod.RisnbrIter = 0; + } + } else { + // The valve is open, check if it should close. + // CLOSING CONDITION: Check existence of a backflow + if (meanFl < 0.) { + eq[cEq].bc[iBc].clsFlgRis = 1; + if (!cm.seq()) { + cm.bcast(cm_mod, &eq[cEq].bc[iBc].clsFlgRis); + } + std::cout << "!!! -- Going from open to close " << std::endl; + com_mod.RisnbrIter = 0; + } + } + + // Check for the status + // If the valve is closed, check the pressure difference, + // if the pressure difference is negative the valve should be open + // -> the status is then not admissible + if (eq[cEq].bc[iBc].clsFlgRis == 1) { + if (eq[cEq].bc[iBc].g < meanP) { + std::cout << "** Not admissible, should be open **" << std::endl; + } + } else { + // If the valve is open, chech the flow, + // if the flow is negative the valve should be closed + // -> the status is then not admissible + if (meanFl < 0.) { + std::cout << "** Not admissible, should be closed **" << std::endl; + } + } + + } + + if (!cm.seq()) { + cm.bcast(cm_mod, &com_mod.RisnbrIter); + } + +} + +} \ No newline at end of file diff --git a/Code/Source/solver/ris.h b/Code/Source/solver/ris.h new file mode 100644 index 000000000..48816fb52 --- /dev/null +++ b/Code/Source/solver/ris.h @@ -0,0 +1,62 @@ +/* Copyright (c) Stanford University, The Regents of the University of California, and others. + * + * All Rights Reserved. + * + * See Copyright-SimVascular.txt for additional details. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject + * to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS + * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER + * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef RIS_H +#define RIS_H + +#include "ComMod.h" + +namespace ris { + +void ris_meanq(ComMod& com_mod, CmMod& cm_mod); +void ris_resbc(ComMod& com_mod, const Array& Yg, const Array& Dg); +void setbc_ris(ComMod& com_mod, const bcType& lBc, const mshType& lM, const faceType& lFa, + const Array& Yg, const Array& Dg); + +void ris_updater(ComMod& com_mod, CmMod& cm_mod); +void ris_status(ComMod& com_mod, CmMod& cm_mod); + +void doassem_ris(ComMod& com_mod, const int d, const Vector& eqN, + const Array3& lK, const Array& lR); + +void doassem_velris(ComMod& com_mod, const int d, const Array& eqN, + const Array3& lK, const Array& lR); + +void clean_r_ris(ComMod& com_mod); +void setbcdir_ris(ComMod& com_mod, Array& lA, Array& lY, Array& lD); + +// TODO: RIS 0D code +void ris0d_bc(ComMod& com_mod, CmMod& cm_mod, const Array& Yg, const Array& Dg); +void ris0d_status(ComMod& com_mod, CmMod& cm_mod); //, const Array& Yg, const Array& Dg); + +}; + +#endif + diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 731b86779..266b99608 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1344,6 +1344,8 @@ void set_bc_neu(ComMod& com_mod, const CmMod& cm_mod, const Array& Yg, c dmsg << "----- iBc " << iBc+1; #endif + if (utils::btest(bc.bType, iBC_Ris0D)) {continue;} + if (utils::btest(bc.bType, iBC_Neu)) { #ifdef debug_set_bc_neu dmsg << "iM: " << iM+1; diff --git a/Code/Source/solver/uris.cpp b/Code/Source/solver/uris.cpp new file mode 100644 index 000000000..bfb539dea --- /dev/null +++ b/Code/Source/solver/uris.cpp @@ -0,0 +1,1242 @@ +/* Copyright (c) Stanford University, The Regents of the University of California, and others. + * + * All Rights Reserved. + * + * See Copyright-SimVascular.txt for additional details. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject + * to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS + * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER + * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include "uris.h" + +#include "all_fun.h" +#include "lhsa.h" +#include "mat_fun.h" +#include "nn.h" +#include "utils.h" +#include "set_bc.h" + +#include "load_msh.h" +#include "vtk_xml.h" +#include "read_msh.h" +#include "VtkData.h" + +namespace uris { + +/// @brief This subroutine computes the mean pressure and flux on the +/// immersed surface +void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris) { + #define n_debug_uris_meanp + #ifdef debug_uris_meanp + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "urisFlag: " << com_mod.urisFlag; + #endif + + using namespace consts; + + auto& cm = com_mod.cm; + auto& eq = com_mod.eq; + auto& uris = com_mod.uris; + auto& msh = com_mod.msh; + auto& uris_obj = uris[iUris]; + + const int nsd = com_mod.nsd; + // const int cEq = com_mod.cEq; + + // auto& An = com_mod.An; + // auto& Ad = com_mod.Ad; + // auto& Dn = com_mod.Dn; + auto& Yn = com_mod.Yn; + + // Let's conpute the mean pressure in the two regions of the fluid mesh + // For the moment let's define a flag IdSubDmn(size the number of elements) + + // Now we can compute the pressure mean on each subdomain + // We need to have a sdf array for each mesh + double Deps = uris_obj.sdf_deps * 2.5; + double volU = 0.0; + double volD = 0.0; + + // Let's compute left side + Array sUPS(1,com_mod.tnNo); + // std::cout << "com_mod.tnNo: " << com_mod.tnNo << std::endl; + // std::cout << "uris_obj.sdf size: " << uris_obj.sdf.size() << std::endl; + + sUPS = 0.0; + for (size_t j = 0; j < sUPS.size(); j++) { + if (uris_obj.sdf(j) >= 0.0 && uris_obj.sdf(j) <= Deps) { + // Reverse the sdf distance for aortic valve + // [HZ] Adjust this value to be more flexible about the box + // if (uris_obj.sdf(j) < 0.0 && uris_obj.sdf(j) >= -Deps) { + sUPS(0,j) = 1.0; + } + } + + for (int iM = 0; iM < com_mod.nMsh; iM++) { + volU += all_fun::integ(com_mod, cm_mod, iM, sUPS); + } + + + // Let's compute right side + Array sDST(1,com_mod.tnNo); + sDST = 0.0; + for (size_t j = 0; j < sDST.size(); j++) { + if (uris_obj.sdf(j) < 0.0 && uris_obj.sdf(j) >= -Deps) { + // Reverse the sdf distance for aortic valve + // if (uris_obj.sdf(j) >= 0.0 && uris_obj.sdf(j) <= Deps) { + sDST(0,j) = 1.0; + } + } + + for (int iM = 0; iM < com_mod.nMsh; iM++) { + volD += all_fun::integ(com_mod, cm_mod, iM, sDST); + } + + // Print volume messages. + if (cm.mas(cm_mod)) { + std::cout << "volume upstream " << volU << " for: " << uris_obj.name << std::endl; + std::cout << "volume downstream " << volD << " for: " << uris_obj.name << std::endl; + } + + double meanPU = 0.0; + double meanPD = 0.0; + + int iEq = 0; + int m = 1; + int s = eq[iEq].s + nsd; + // int e = s + m - 1; + + Array tmpV(maxNSD, com_mod.tnNo); + for (int j = 0; j < Yn.ncols(); j++) { + tmpV(0,j) = Yn(s,j)*sUPS(j); + } + for (int iM = 0; iM < com_mod.nMsh; iM++) { + meanPU += all_fun::integ(com_mod, cm_mod, iM, tmpV); + } + meanPU = meanPU / volU; + + for (int j = 0; j < Yn.ncols(); j++) { + tmpV(0,j) = Yn(s,j)*sDST(j); + } + for (int iM = 0; iM < com_mod.nMsh; iM++) { + meanPD += all_fun::integ(com_mod, cm_mod,iM, tmpV); + } + meanPD = meanPD / volD; + + uris_obj.meanPU = uris_obj.relax_factor * meanPU + + (1.0 - uris_obj.relax_factor) * uris_obj.meanPU; + uris_obj.meanPD = uris_obj.relax_factor * meanPD + + (1.0 - uris_obj.relax_factor) * uris_obj.meanPD; + + if (cm.mas(cm_mod)) { + std::cout << "mean P upstream " << meanPU << " " << uris_obj.meanPU + << " for: " << uris_obj.name << std::endl; + std::cout << "mean P downstream " << meanPD << " " << uris_obj.meanPD + << " for: " << uris_obj.name << std::endl; + } + + // If the uris has passed the closing state + if (uris_obj.cnt > uris_obj.DxClose.nrows()) { + if (uris_obj.meanPD > uris_obj.meanPU) { + uris_obj.cnt = 1; + uris_obj.clsFlg = false; + com_mod.urisActFlag = true; + if (cm.mas(cm_mod)) { + std::cout << "** Set urisCloseFlag to FALSE for: " + << uris_obj.name << std::endl; + } + } + } + if (cm.mas(cm_mod)) { + std::cout << "urisCloseFlag is: " << uris_obj.clsFlg << " for: " + << uris_obj.name << std::endl; + } + +} + +/// @brief This subroutine computes the mean velocity in the fluid elements +/// near the immersed surface +void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris) { + #define n_debug_uris_meanv + #ifdef debug_uris_meanv + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "urisFlag: " << com_mod.urisFlag; + #endif + + using namespace consts; + + auto& cm = com_mod.cm; + auto& eq = com_mod.eq; + auto& uris = com_mod.uris; + auto& msh = com_mod.msh; + auto& uris_obj = uris[iUris]; + + const int nsd = com_mod.nsd; + // const int cEq = com_mod.cEq; + + // auto& An = com_mod.An; + // auto& Ad = com_mod.Ad; + // auto& Dn = com_mod.Dn; + auto& Yn = com_mod.Yn; + + // Let's compute the neighboring region below the valve normal. When + // the valve is open, this region should roughly be valve oriface. + int iEq = 0; + + double Deps = uris_obj.sdf_deps; + Array sImm(1, com_mod.tnNo); + sImm = 0.0; + double volI = 0.0; + + for (int i = 0; i < com_mod.tnNo; i++) { + if (uris_obj.sdf(i) <= -Deps) { + // Reverse the sdf distance for aortic valve + // if (uris_obj.sdf(i) >= Deps) { + sImm(0,i) = 1.0; + } + } + + for (int iM = 0; iM < com_mod.nMsh; iM++) { + volI += all_fun::integ(com_mod, cm_mod, iM, sImm); + } + if (cm.mas(cm_mod)) { + std::cout << "volume inside " << volI << " for: " << uris_obj.name << std::endl; + } + + int m = nsd; + int s = eq[iEq].s; + // int e = s + m - 1; + + Array tmpV(maxNSD, com_mod.tnNo); + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < Yn.ncols(); j++) { + tmpV(i,j) = Yn(s+i,j)*sImm(0,j); + } + } + + Array tmpVNrm(1, com_mod.tnNo); + for (int i = 0; i < com_mod.tnNo; i++) { + tmpVNrm(0,i) = tmpV(0,i) * uris_obj.nrm(0) + + tmpV(1,i) * uris_obj.nrm(1) + + tmpV(2,i) * uris_obj.nrm(2); + } + + double meanV = 0.0; + for (int iM = 0; iM < com_mod.nMsh; iM++) { + meanV += all_fun::integ(com_mod, cm_mod, iM, tmpVNrm)/volI; + } + + if (cm.mas(cm_mod)) { + std::cout << "mean velocity: " << meanV << " for: " << uris_obj.name << std::endl; + } + + // If the uris has passed the open state + if (uris_obj.cnt > uris_obj.DxOpen.nrows()) { + if (meanV < 0.0) { + uris_obj.cnt = 1; + uris_obj.clsFlg = true; + com_mod.urisActFlag = true; + if (cm.mas(cm_mod)) { + std::cout << "** Set urisCloseFlag to TRUE for: " + << uris_obj.name << std::endl; + } + } + } + if (cm.mas(cm_mod)) { + std::cout << "urisCloseFlag is: " << uris_obj.clsFlg << " for: " + << uris_obj.name << std::endl; + } + +} + +/// @brief This subroutine computes the displacement of the immersed +/// surface with fem projection +void uris_update_disp(ComMod& com_mod, CmMod& cm_mod) { + #define n_debug_uris_update_disp + #ifdef debug_uris_update_disp + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + // using namespace consts; + + auto& cm = com_mod.cm; + // auto& eq = com_mod.eq; + auto& uris = com_mod.uris; + auto& msh = com_mod.msh; + int nUris = com_mod.nUris; + + const int nsd = com_mod.nsd; + + // For each point in the immersed surface we need to localize it + // = find the fluid element that contains the node + // Since the fluid element could be on another processor, we need to + // gather the displacement values at the end + // [FK] it's probably better to save the element ids so that we don't + // have to run the search every time step, only during open or close + + Array localYd, xl, Nxi; + Vector N; + Vector xp(nsd), xi(nsd), d(nsd); + bool fl; + + for (int iUris = 0; iUris < nUris; iUris++) { + auto& uris_obj = uris[iUris]; + uris_find_tetra(com_mod, cm_mod, iUris); + localYd.resize(nsd, uris_obj.tnNo); + localYd = 0.0; + for (int nd = 0; nd < uris_obj.tnNo; nd++) { + int jM = uris_obj.elemId(0, nd); + auto& mesh = msh[jM]; + // If the fluid mesh element is not on the current proc + if (jM == -1) {continue;} + + int iEln = uris_obj.elemId(1, nd); + Vector xp = uris_obj.x.col(nd); + xl.resize(nsd, mesh.eNoN); + N.resize(mesh.eNoN); + Nxi.resize(nsd, mesh.eNoN); + + for (int a = 0; a < mesh.eNoN; a++) { + int Ac = mesh.IEN(a, iEln); + for (int i = 0; i < nsd; i++) { + xl(i, a) = com_mod.x(i, Ac); + } + } + // Get displacement + // Localize p inside the parent element + nn::get_xi(nsd, mesh.eType, mesh.eNoN, xl, xp, xi, fl); + if (!fl) { + if (cm.mas(cm_mod)) { + std::cout << "[WARNING] URIS get_xi not converging!" << std::endl; + } + } + // evaluate N at xi + nn::get_gnn(nsd, mesh.eType, mesh.eNoN, xi, N, Nxi); + // use this to compute disp al node xp + d = 0.0; + for (int a = 0; a < mesh.eNoN; a++) { + int Ac = mesh.IEN(a, iEln); + //We have to use Do because Dn contains the result coming from the solid + d(0) += N(a)*com_mod.Do(nsd+1, Ac); + d(1) += N(a)*com_mod.Do(nsd+2, Ac); + d(2) += N(a)*com_mod.Do(nsd+3, Ac); + } + // update uris disp + localYd.set_col(nd, d); + } + MPI_Allreduce(localYd.data(), uris_obj.Yd.data(), uris_obj.tnNo*nsd, + cm_mod::mpreal, MPI_SUM, cm.com()); + + for (int nd = 0; nd < uris_obj.tnNo; nd++) { + double divisor = static_cast(std::max(1, uris_obj.elemCounter(nd))); + // Vector Yd_vec = uris_obj.Yd.col(nd) / div; + // uris_obj.Yd.set_col(nd, Yd_vec); + for (int i = 0; i < nsd; i++) { + uris_obj.Yd(i,nd) /= divisor; + } + } + + } + +} + +/// @brief This subroutine computes the tetrahedral elements +void uris_find_tetra(ComMod& com_mod, CmMod& cm_mod, const int iUris) { + #define n_debug_uris_find_tetra + #ifdef debug_uris_find_tetra + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + auto& cm = com_mod.cm; + auto& uris = com_mod.uris; + auto& uris_obj = uris[iUris]; + const int nsd = com_mod.nsd; + + // We need to check if the valve needs to move + int cnt; + if (!uris_obj.clsFlg) { + cnt = std::min(uris_obj.cnt, uris_obj.DxOpen.nrows()); + } else { + cnt = std::min(uris_obj.cnt, uris_obj.DxClose.nrows()); + } + + if (uris_obj.elemId.allocated() && cnt < uris_obj.cnt) { + return; + } + + // For each point in the immersed surface we need to localize it + // = find the fluid element that contains the node + // Since the fluid element could be on another processor, we need to + // gather the displacement values at the end + // [FK] it's probably better to save the element ids so that we don't + // have to run the search every time step, only during open or close + + bool ultra = true; + if (!uris_obj.elemId.allocated()) { + uris_obj.elemId.resize(2, uris_obj.tnNo); + } + if (!uris_obj.elemCounter.allocated()) { + uris_obj.elemCounter.resize(uris_obj.tnNo); + } + Vector local_counter(uris_obj.tnNo); + local_counter = 0; + uris_obj.elemId = -1; + uris_obj.elemCounter = 0; + int flag; + Array xl; + for (int nd = 0; nd < uris_obj.tnNo; nd++) { + flag = 0; + // Check if we were able to find the tetra. + // [FK] if not, the tetra is on another processor + Vector xp = uris_obj.x.col(nd); + bool found = false; + + for (int jM = 0; jM < com_mod.nMsh && !found; jM++) { + auto& mesh = com_mod.msh[jM]; + xl.resize(nsd, mesh.eNoN); + for (int iEln = 0; iEln < mesh.nEl && !found; iEln++) { + for (int a = 0; a < mesh.eNoN; a++) { + int Ac = mesh.IEN(a, iEln); + for(int i = 0; i < nsd; i++) { + xl(i,a) = com_mod.x(i,Ac); + } + } + inside_tet(com_mod, mesh.eNoN, xp, xl, flag, ultra); + if (flag == 1) { + uris_obj.elemId(0, nd) = jM; + uris_obj.elemId(1, nd) = iEln; + local_counter(nd) += 1; + found = true; + } + } + } + } + + MPI_Allreduce(local_counter.data(), uris_obj.elemCounter.data(), + uris_obj.tnNo, cm_mod::mpint, MPI_SUM, cm.com()); + +} + + +/// @brief This subroutine check if a node is inside a tetrahedron +void inside_tet(ComMod& com_mod, int& eNoN, Vector& xp, + Array& xl, int& flag, bool ext) { + #define n_debug_inside_tet + #ifdef debug_inside_tet + DebugMsg dmsg(__func__, com_mod.cm.idcm()); + dmsg.banner(); + #endif + + auto& uris = com_mod.uris; + const int nsd = com_mod.nsd; + Vector minb(nsd); + Vector maxb(nsd); + + // Create a bounding box around of the current solid location + // [FK]: Hard coded BBox?? This is going to cause problem if scale changes + for (int i = 0; i < nsd; i++) { + double min_val = std::numeric_limits::max(); + double max_val = std::numeric_limits::lowest(); + for (int j = 0; j < eNoN; j++) { + double val_minus = xl(i,j) - 0.1; + double val_plus = xl(i,j) + 0.1; + if (val_minus < min_val) {min_val = val_minus;} + if (val_plus > max_val) {max_val = val_plus;} + } + minb(i) = min_val; + maxb(i) = max_val; + } + + // Is the node inside the BBox? + bool inside = true; + for (int i = 0; i < nsd; ++i) { + if (xp(i) < minb(i) || xp(i) > maxb(i)) { + inside = false; + break; + } + } + + flag = 0; + if (inside) { + flag = in_poly(xp, xl, ext); + } +} + + +/// @brief Read the URIS mesh separately +void uris_read_msh(Simulation* simulation) { + #define n_debug_uris_read_msh + #ifdef debug_uris_read_msh + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); + dmsg.banner(); + #endif + + auto& com_mod = simulation->get_com_mod(); + auto& uris = com_mod.uris; + const int nsd = com_mod.nsd; + + com_mod.urisFlag = true; + com_mod.urisActFlag = true; + + auto param = simulation->parameters.URIS_mesh_parameters[0]; + com_mod.urisRes = param->resistance(); + com_mod.urisResClose = param->resistance_close(); + + std::cout << "URIS resistance: " << com_mod.urisRes << std::endl; + std::cout << "URIS resistance when the valve is closed: " << com_mod.urisResClose << std::endl; + + int nUris = simulation->parameters.URIS_mesh_parameters.size(); + com_mod.nUris = nUris; + + std::cout << "Number of immersed surfaces for uris: " << nUris << std::endl; + uris.resize(nUris); + + for (int iUris = 0; iUris < nUris; iUris++) { + auto param = simulation->parameters.URIS_mesh_parameters[iUris]; + auto& uris_obj = uris[iUris]; + uris_obj.name = param->name(); + std::cout << "** Reading URIS mesh: " << uris_obj.name << std::endl; + + uris_obj.scF = param->mesh_scale_factor(); + uris_obj.nFa = param->URIS_face_parameters.size(); + uris_obj.msh.resize(uris_obj.nFa); + uris_obj.nrm.resize(nsd); + Array gX(0,0); + + std::string positive_flow_normal_file_path = param->positive_flow_normal_file_path(); + // [HZ] Need to read flow normal file (*.dat) into uris_obj.nrm + // lPtr => lPM%get(fTmp, "Positive flow normal file") + // fid = fTmp%open() + // READ (fid,*) uris(iUris)%nrm(:) + // CLOSE (fid) + std::ifstream file_stream; + file_stream.open(positive_flow_normal_file_path); + if (!file_stream.is_open()) { + throw std::runtime_error("Failed to open the open positive flow normal file '" + + positive_flow_normal_file_path + "'."); + } + for (int i = 0; i < nsd; i++) { + file_stream >> uris_obj.nrm(i); + } + file_stream.close(); + + uris_obj.sdf_deps = param->thickness(); + uris_obj.sdf_deps_close = param->close_thickness(); + uris_obj.clsFlg = param->valve_starts_as_closed(); + + // uris_obj.tnNo = 0; + for (int iM = 0; iM < uris_obj.nFa; iM++) { + // Set as shell + auto mesh_param = param->URIS_face_parameters[iM]; + auto& mesh = uris_obj.msh[iM]; + mesh.lShl = true; + mesh.name = mesh_param->name(); + std::cout << "-- Reading URIS face: " << mesh.name << std::endl; + + // Read mesh nodal coordinates and element connectivity. + uris_read_sv(simulation, mesh, mesh_param); + + // // [HZ] eType is not NA, neglecting this for now + // IF (uris(iUris)%msh(iM)%eType .EQ. eType_NA) THEN + // CALL READCCNE(lPN, uris(iUris)%msh(iM)) + // END IF + // IF (uris(iUris)%msh(iM)%eType .EQ. eType_NA) THEN + // CALL READNRB(lPN, uris(iUris)%msh(iM)) + // END IF + // IF (uris(iUris)%msh(iM)%eType .EQ. eType_NA) THEN + // CALL READGAMBIT(lPN, uris(iUris)%msh(iM)) + // END IF + // IF (uris(iUris)%msh(iM)%eType .EQ. eType_NA) THEN + // err = " Failed to identify format of the uris mesh" + // END IF + + std::cout << "Number of uris nodes: " << mesh.gnNo << std::endl; + std::cout << "Number of uris elements: " << mesh.gnEl << std::endl; + + // Read valve motion: note that this motion is defined on the + // reference configuration + std::string open_motion_file_path = mesh_param->open_motion_file_path(); + // std::ifstream file_stream; + file_stream.open(open_motion_file_path); + if (!file_stream.is_open()) { + throw std::runtime_error("Failed to open the open motion file '" + + open_motion_file_path + "'."); + } + int dispNtOpen, dispNnOpen; + file_stream >> dispNtOpen >> dispNnOpen; + // std::cout << "dispNtOpen: " << dispNtOpen << std::endl; + // std::cout << "dispNnOpen: " << dispNnOpen << std::endl; + + if (dispNnOpen != mesh.gnNo) { + throw std::runtime_error("Mismatch in node numbers between URIS mesh and displacements."); + } + + Array3 dispOpen(dispNtOpen, nsd, dispNnOpen); + for (int t = 0; t < dispNtOpen; t++) { + for (int a = 0; a < dispNnOpen; a++) { + for (int i = 0; i < nsd; i++) { + file_stream >> dispOpen(t,i,a); + // std::cout << "dispOpen: " << dispOpen(t,n,a) << std::endl; + } + } + } + file_stream.close(); + + std::string close_motion_file_path = mesh_param->close_motion_file_path(); + // std::ifstream file_stream; + file_stream.open(close_motion_file_path); + if (!file_stream.is_open()) { + throw std::runtime_error("Failed to open the close motion file '" + + close_motion_file_path + "'."); + } + int dispNtClose, dispNnClose; + file_stream >> dispNtClose >> dispNnClose; + // std::cout << "dispNtClose: " << dispNtClose << std::endl; + // std::cout << "dispNnClose: " << dispNnClose << std::endl; + + if (dispNnClose != mesh.gnNo) { + throw std::runtime_error("Mismatch in node numbers between URIS mesh and displacements."); + } + + Array3 dispClose(dispNtClose, nsd, dispNnClose); + for (int t = 0; t < dispNtClose; t++) { + for (int a = 0; a < dispNnClose; a++) { + for (int i = 0; i < nsd; i++) { + file_stream >> dispClose(t,i,a); + // std::cout << "dispClose: " << dispClose(t,n,a) << std::endl; + } + } + } + file_stream.close(); + + // To scale the mesh, while attaching x to gX + int a = uris_obj.tnNo + mesh.gnNo; + // std::cout << "uris obj tnNo: " << uris_obj.tnNo << std::endl; + // std::cout << "mesh gnNo: " << mesh.gnNo << std::endl; + // std::cout << "mesh x size: " << mesh.x.nrows() << ", " << mesh.x.ncols() << std::endl; + + if (iM == 0) { + gX.resize(nsd, a); + gX = 0.0; + uris_obj.DxOpen.resize(dispNtOpen, nsd, a); + uris_obj.DxClose.resize(dispNtClose, nsd, a); + } else{ + Array tmpX(nsd, uris_obj.tnNo); + tmpX = gX; + gX.resize(nsd, a); + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < uris_obj.tnNo; j++) { + gX(i,j) = tmpX(i,j); + } + } + // Move data for open + Array3 tmpDxOpen(dispNtOpen, nsd, uris_obj.tnNo); + tmpDxOpen = uris_obj.DxOpen; + uris_obj.DxOpen.resize(dispNtOpen, nsd, a); + for (int k = 0; k < dispNtOpen; k++) { + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < uris_obj.tnNo; j++) { + uris_obj.DxOpen(k,i,j) = tmpDxOpen(k,i,j); + } + } + } + // Move data for open + Array3 tmpDxClose(dispNtClose, nsd, uris_obj.tnNo); + tmpDxClose = uris_obj.DxClose; + uris_obj.DxClose.resize(dispNtClose, nsd, a); + for (int k = 0; k < dispNtClose; k++) { + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < uris_obj.tnNo; j++) { + uris_obj.DxClose(k,i,j) = tmpDxClose(k,i,j); + } + } + } + } + + for (int i = 0; i < nsd; i++) { + for (int j = uris_obj.tnNo; j < a; j++) { + gX(i,j) = mesh.x(i,j-uris_obj.tnNo) * uris_obj.scF; + } + } + for (int k = 0; k < dispNtOpen; k++) { + for (int i = 0; i < nsd; i++) { + for (int j = uris_obj.tnNo; j < a; j++) { + uris_obj.DxOpen(k,i,j) = dispOpen(k,i,j-uris_obj.tnNo); + } + } + } + for (int k = 0; k < dispNtClose; k++) { + for (int i = 0; i < nsd; i++) { + for (int j = uris_obj.tnNo; j < a; j++) { + uris_obj.DxClose(k,i,j) = dispClose(k,i,j-uris_obj.tnNo); + } + } + } + uris_obj.tnNo = a; + // mesh.x.clear(); + // dispOpen.clear(); + // dispClose.clear(); + } + uris_obj.x.resize(nsd, uris_obj.tnNo); + uris_obj.x = gX; + uris_obj.Yd.resize(nsd, uris_obj.tnNo); + uris_obj.Yd = 0.0; + // gX.clear(); + + // Setting mesh.gN, mesh.lN parameter + int b = 0; + for (int iM = 0; iM < uris_obj.nFa; iM++) { + auto& mesh = uris_obj.msh[iM]; + mesh.nNo = mesh.gnNo; + mesh.gN.resize(mesh.nNo); + mesh.gN = 0; + mesh.lN.resize(uris_obj.tnNo); + mesh.lN = 0; + for (int a = 0; a < mesh.nNo; a++) { + mesh.gN(a) = b; + mesh.lN(b) = a; + b++; + } + } + if (b != uris_obj.tnNo) { + throw std::runtime_error("Mismatch in uris.tnNo. Correction needed."); + } + + // Remap msh%gIEN array + for (int iM = 0; iM < uris_obj.nFa; iM++) { + auto& mesh = uris_obj.msh[iM]; + mesh.nEl = mesh.gnEl; + mesh.IEN.resize(mesh.eNoN, mesh.nEl); + for (int e = 0; e < mesh.nEl; e++) { + for (int a = 0; a < mesh.eNoN; a++) { + int Ac = mesh.gIEN(a,e); + Ac = mesh.gN(Ac); + mesh.IEN(a,e) = Ac; + } + } + mesh.gIEN.clear(); + } + + if (uris_obj.nFa > 0) { + std::string msg = "Total number of uris nodes: " + std::to_string(uris_obj.tnNo); + std::cout << msg << std::endl; + int total_nel = 0; + for (int iM = 0; iM < uris_obj.nFa; iM++) { + total_nel += uris_obj.msh[iM].nEl; + } + msg = "Total number of uris elements: " + std::to_string(total_nel); + std::cout << msg << std::endl; + } + } + std::cout << "URIS mesh data imported successfully." << std::endl; + +} + +/// @brief Write URIS solution to a vtu file +void uris_write_vtus(ComMod& com_mod) { + #define n_debug_uris_write_vtus + #ifdef debug_uris_write_vtus + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); + dmsg.banner(); + #endif + + using namespace consts; + + auto& uris = com_mod.uris; + const int nsd = com_mod.nsd; + const int nUris = com_mod.nUris; + + // we plot coord + displacement + int nOut = 2; + int outDof = nOut * nsd; + std::vector outNames(nOut); + Vector outS(nOut+1); + + // Prepare all solultions in to dataType d + std::vector d; + for (int iUris = 0; iUris < nUris; iUris++) { + auto& uris_obj = uris[iUris]; + d.resize(uris_obj.nFa); + int nNo = 0; + int nEl = 0; + for (int iM = 0; iM < uris_obj.nFa; iM++) { + auto& mesh = uris_obj.msh[iM]; + int cOut = 0; + outS(cOut) = 0; // [HZ] Need to check this if it's 1 or 0 + outS(cOut+1) = nsd; + // outNames[cOut] = ""; + + // outS = [0, 3] + + if (mesh.eType == ElementType::NRB) { + throw std::runtime_error("Outputs for NURBS data is under development."); + } + d[iM].nNo = mesh.nNo; + d[iM].nEl = mesh.nEl; + d[iM].eNoN = mesh.eNoN; + d[iM].vtkType = mesh.vtkType; + d[iM].x.resize(outDof, mesh.nNo); + d[iM].IEN.resize(mesh.eNoN, mesh.nEl); + + for (int a = 0; a < mesh.nNo; a++) { + int Ac = mesh.gN(a); + for (int i = 0; i < nsd; i++) { + d[iM].x(i,a) = uris_obj.x(i,Ac); + } + } + for (int e = 0; e < mesh.nEl; e++) { + for (int i = 0; i < mesh.eNoN; i++) { + d[iM].IEN(i,e) = mesh.IEN(i,e); + } + } + + int l = nsd; // l = 3 + int s = 0; + int e = s + l; // e = 3 + + cOut += 1; // cOut = 1 + int is = outS(cOut); // is = outS[1] = 3 + int ie = is + l; // ie = 3 + 3 = 6 + outS(cOut+1) = ie; // outS[2] = 7, outS = [0, 3, 6] + outNames[0] = "coordinates"; + outNames[1] = "URIS_displacement"; + + for (int a = 0; a < mesh.nNo; a++) { + int Ac = mesh.gN(a); + for (int i = 0; i < nsd; i++) { + d[iM].x(is+i,a) = uris_obj.Yd(s+i,Ac); // [HZ] Need to check this + } + } + + nNo += mesh.nNo; + nEl += mesh.nEl; + } + + // Writing to vtu file (master only) + char fName_num[100]; + if (com_mod.cTS >= 1000) { + sprintf(fName_num, "%d", com_mod.cTS); + } else { + sprintf(fName_num, "%03d", com_mod.cTS); + } + std::string fName = com_mod.saveName + "_uris_" + uris_obj.name + "_" + fName_num + ".vtu"; + + auto vtk_writer = VtkData::create_writer(fName); + // Writing the position data + int iOut = 0; + int s = outS(iOut); + // int e = outS(iOut+1) - 1; + int nSh = 0; + Array tmpV(maxNSD, nNo); + tmpV = 0.0; + for (int iM = 0; iM < uris_obj.nFa; iM++) { + for (int a = 0; a < d[iM].nNo; a++) { + for (int i = 0; i < nsd; i++) { + tmpV(i,a+nSh) = d[iM].x(s+i,a); + } + } + nSh += d[iM].nNo; + } + vtk_writer->set_points(tmpV); + + // Writing the connectivity data + Array tmpI; + for (int iM = 0; iM < uris_obj.nFa; iM++) { + tmpI.resize(d[iM].eNoN, d[iM].nEl); + for (int e = 0; e < d[iM].nEl; e++) { + for (int i = 0; i < d[iM].eNoN; i++) { + tmpI(i,e) = d[iM].IEN(i,e); + } + } + vtk_writer->set_connectivity(nsd, tmpI); + } + + // Writing all solutions + for (int iOut = 0; iOut < nOut; iOut++) { + int s = outS(iOut); + int e = outS(iOut+1); + int l = e - s; + tmpV.resize(l, nNo); + tmpV = 0.0; + int nSh = 0; + for (int iM = 0; iM < uris_obj.nFa; iM++) { + for (int a = 0; a < d[iM].nNo; a++) { + for (int i = 0; i < nsd; i++) { + tmpV(i,a+nSh) = d[iM].x(s+i,a); + } + } + nSh += d[iM].nNo; + } + vtk_writer->set_point_data(outNames[iOut], tmpV); + } + + vtk_writer->write(); + + delete vtk_writer; + + } + +} + +/// @brief Checks if a probe lies inside or outside an immersed boundary +void uris_calc_sdf(ComMod& com_mod) { + #define n_debug_uris_calc_sdf + #ifdef debug_uris_calc_sdf + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); + dmsg.banner(); + #endif + + auto& cm = com_mod.cm; + auto& uris = com_mod.uris; + const int nsd = com_mod.nsd; + const int nUris = com_mod.nUris; + + Array xXi(nsd, nsd-1); + + for (int iUris = 0; iUris < nUris; iUris++) { + // We need to check if the valve needs to move + auto& uris_obj = uris[iUris]; + int cnt = 0; + if (!uris_obj.clsFlg) { + cnt = std::min(uris_obj.cnt, uris_obj.DxOpen.nrows()); + for (int i = 0; i < uris_obj.x.nrows(); i++) { + for (int j = 0; j < uris_obj.x.ncols(); j++) { + uris_obj.x(i,j) = uris_obj.DxOpen(cnt-1,i,j); + } + } + } else { + cnt = std::min(uris_obj.cnt, uris_obj.DxClose.nrows()); + for (int i = 0; i < uris_obj.x.nrows(); i++) { + for (int j = 0; j < uris_obj.x.ncols(); j++) { + uris_obj.x(i,j) = uris_obj.DxClose(cnt-1,i,j); + } + } + } + + // if (uris_obj.sdf.allocated() && cnt < uris_obj.cnt) {continue;} + if (uris_obj.sdf.size() > 0 && cnt < uris_obj.cnt) {continue;} + + int max_eNoN = 0; + for (int iM = 0; iM < uris_obj.nFa; iM++) { + auto& mesh = uris_obj.msh[iM]; + if (mesh.eNoN > max_eNoN) { + max_eNoN = mesh.eNoN; + } + } + + Array lX(nsd, uris_obj.msh[1].eNoN); + // if (!uris_obj.sdf.allocated()) { + if (uris_obj.sdf.size() <= 0) { + uris_obj.sdf.resize(com_mod.tnNo); + uris_obj.sdf = 0.0; + } + + if (cm.idcm() == 0) { + std::cout << "Recomputing SDF for " << uris_obj.name << std::endl; + } + uris_obj.sdf = uris_obj.sdf_default; + + // Each time when the URIS moves (open/close), we need to + // recompute the signed distance function. + // Find the bounding box of the valve, the BBox will be 10% larger + // than the actual valve. + Vector minb(nsd); + Vector maxb(nsd); + Vector extra(nsd); + for (int i = 0; i < nsd; i++) { + minb(i) = std::numeric_limits::max(); + maxb(i) = std::numeric_limits::lowest(); + } + + // For each coordinate dimension, find the minimum and maximum in uris_obj.x. + double extra_val = 0.1; // [HZ] The BBox is 10% larger than the actual valve, default is 0.1 + for (int i = 0; i < nsd; i++) { + for (int j = 0; j < uris_obj.x.ncols(); j++) { + double val = uris_obj.x(i,j); + if (val < minb(i)) + minb(i) = val; + if (val > maxb(i)) + maxb(i) = val; + } + extra(i) = (maxb(i) - minb(i)) * extra_val; + } + + // The SDF is computed on the reference configuration, which + // means that the valves will be morphed based on the fluid mesh + // motion. If the fluid mesh stretches near the valve, the valve + // leaflets will also be streched. Note that + // this is a simplifying assumption. + Vector xp(nsd); + for (int ca = 0; ca < com_mod.tnNo; ca++) { + double minS = std::numeric_limits::max(); + for (int i = 0; i < nsd; i++) { + xp(i) = com_mod.x(i,ca); + } + // Is the node inside the BBox? + bool inside = true; + for (int i = 0; i < nsd; i++) { + if (xp(i) < (minb(i) - extra(i)) || xp(i) > (maxb(i) + extra(i))) { + inside = false; + break; + } + } + if (inside) { + // This point is inside the BBox + // Find the closest URIS face centroid + int Ec = -1; + int jM = -1; + Vector xb(nsd); + for (int iM = 0; iM < uris_obj.nFa; iM++) { + auto& mesh = uris_obj.msh[iM]; + for (int e = 0; e < mesh.nEl; e++) { + xb = 0.0; + for (int a = 0; a < mesh.eNoN; a++) { + int Ac = mesh.IEN(a,e); + for (int i = 0; i < nsd; i++) { + xb(i) += uris_obj.x(i,Ac); + } + } + for (int i = 0; i < nsd; i++) { + xb(i) /= static_cast(mesh.eNoN); + } + double dS = 0.0; + for (int i = 0; i < nsd; i++) { + dS += (xp[i] - xb[i]) * (xp[i] - xb[i]); + } + dS = std::sqrt(dS); + + if (dS < minS) { + minS = dS; + Ec = e; + jM = iM; + } + } + } + + // We also need to compute the sign (above or below the valve). + // Compute the element normal + auto& mesh = uris_obj.msh[jM]; + xXi = 0.0; + lX = 0.0; + xb = 0.0; + for (int a = 0; a < mesh.eNoN; a++) { + int Ac = mesh.IEN(a,Ec); + for (int i = 0; i < nsd; i++) { + xb(i) += uris_obj.x(i,Ac); + lX(i,a) = uris_obj.x(i,Ac); + } + } + for (int i = 0; i < nsd; i++) { + xb(i) /= static_cast(mesh.eNoN); + } + + for (int a = 0; a < mesh.eNoN; a++) { + for (int i = 0; i < nsd - 1; i++) { + double factor = mesh.Nx(i,a,0); + for (int j = 0; j < nsd; j++) + xXi(j,i) += lX(j,a) * factor; + } + } + + auto nV = utils::cross(xXi); + auto Jac = sqrt(utils::norm(nV)); + nV = nV / Jac; + auto dotP = utils::norm(xp-xb, nV); + + // if (dotP < 0.0) { + // dotP = -1.0; + // } else { + // dotP = 1.0; + // } + + // [HZ] Improved implementation for SDF sign + if (uris_obj.clsFlg) { + auto dot_nrm = utils::norm(xp-xb, uris_obj.nrm); + if (dot_nrm < 0.0 && dotP < 0.0) { + dotP = -1.0; + } else { + dotP = 1.0; + } + } else { + if (dotP < 0.0) { + dotP = -1.0; + } else { + dotP = 1.0; + } + } + + uris_obj.sdf[ca] = dotP * minS; + } + } + } +} + + +/// @brief Create data for a mesh. +/// +/// Replicates Fortran READSV subroutine defined in LOADMSH.f. +/// +/// SUBROUTINE READSV(list, lM) +// + +void uris_read_sv(Simulation* simulation, mshType& mesh, const URISFaceParameters* mesh_param) { + #define n_dbg_read_sv + #ifdef dbg_read_sv + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "Mesh name: " << mesh_name; + dmsg << "Mesh path: " << mesh_path; + dmsg << "mesh.lShl: " << mesh.lShl; + #endif + + auto mesh_path = mesh_param->face_file_path(); + auto mesh_name = mesh_param->name(); + // Read in volume mesh. + vtk_xml::read_vtu(mesh_path, mesh); + + // Check that the input number of spatial dimensions is consistent + // with the types of elements defined for the simulation mesh. + // + int nsd = simulation->com_mod.nsd; + int elem_dim = consts::element_dimension.at(mesh.eType); + auto elem_type = consts::element_type_to_string.at(mesh.eType); + + if (mesh.lShl) { + if (nsd == 1) { + throw std::runtime_error("The number of spatial dimensions (" + std::to_string(nsd) + + ") is not consistent with the mesh '" + mesh.name + "' which contains shell elements."); + } + + } else if (!mesh.lFib) { + if (elem_dim != nsd) { + throw std::runtime_error("The number of spatial dimensions (" + std::to_string(nsd) + + ") is not consistent with the mesh '" + mesh.name + "' which contains " + elem_type + " elements."); + } + } + + // Set mesh element properites for the input element type. + nn::select_ele(simulation->com_mod, mesh); + + // Check the mesh element node ordering. + // + // Note: This may change element node ordering. + // + auto &com_mod = simulation->get_com_mod(); + if (com_mod.ichckIEN) { + read_msh_ns::check_ien(simulation, mesh); + } +} + + +/// @brief This routine gives the distance between two points + +int in_poly(Vector& P, Array& P1, bool ext) { + int nd = P1.nrows(); + bool flag = true; + int inpoly = 0; + + Vector N(nd); + N = 0.0; + + if (nd == 2) { + for (int i = 0; i <= nd; i++) { + // compute normal in 2D for P2-P1 P3-P1 + if (i != nd) { + N(0) = P1(1,i) - P1(1,i+1); + N(1) = P1(0,i+1) - P1(0,i); + } else { + N(0) = P1(1,i) - P1(1,0); + N(1) = P1(0,0) - P1(0,i); + } + // test dot product between P-P1 and the normals + double dotP = N(0)*(P(0)-P1(0,i)) + N(1)*(P(1)-P1(1,i)); + if (dotP < 0.0) { + flag = false; + break; + } + } + if (flag) { + inpoly = 1; + } + } else { + Vector v1 = P1.col(0); + Vector v2 = P1.col(1); + Vector v3 = P1.col(2); + Vector v4 = P1.col(3); + int s1 = same_side(v1, v2, v3, v4, P, ext); + int s2 = same_side(v2, v3, v4, v1, P, ext); + int s3 = same_side(v3, v4, v1, v2, P, ext); + int s4 = same_side(v4, v1, v2, v3, P, ext); + inpoly = (s1 + s2 + s3 + s4) / 4; + } + + return inpoly; +} + + +/// @brief Chech if a point is on the same side of anotehr point wrt a triangle in 3D +int same_side(Vector& v1, Vector& v2, Vector& v3, + Vector& v4, Vector& p, bool ext) { + #define n_dbg_read_sv + #ifdef dbg_read_sv + DebugMsg dmsg(__func__, simulation->com_mod.cm.idcm()); + dmsg.banner(); + dmsg << "checking same side"; + #endif + + int sameside = 0; + double eps = 2.0e-4; + + Vector v21 = v2 - v1; + Vector v31 = v3 - v1; + Array V(3,2); + V.set_col(0, v21); + V.set_col(1, v31); + Vector v41 = v4 - v1; + Vector vp1 = p - v1; + Vector N = utils::cross(V); + double dotV4 = utils::norm(N, v41); + double dotP = utils::norm(N, vp1); + // check if P and P4 are from the same side + // int sn = utils::sign(dotP); + // if (sn == 1) { + if ((dotP >= 0 && dotV4 >= 0) || (dotP < 0 && dotV4 < 0)) { + sameside = 1; + } + + // If it is not, check if it is on any face, it might on the face + if (ext && sameside != 1) { + if (std::fabs(dotP) <= eps) { + sameside = 1; + } + } + + return sameside; +} + +} \ No newline at end of file diff --git a/Code/Source/solver/uris.h b/Code/Source/solver/uris.h new file mode 100644 index 000000000..cb3b0856d --- /dev/null +++ b/Code/Source/solver/uris.h @@ -0,0 +1,66 @@ +/* Copyright (c) Stanford University, The Regents of the University of California, and others. + * + * All Rights Reserved. + * + * See Copyright-SimVascular.txt for additional details. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject + * to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS + * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER + * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#ifndef URIS_H +#define URIS_H + +#include "ComMod.h" +#include "Simulation.h" + +namespace uris { + +void uris_meanp(ComMod& com_mod, CmMod& cm_mod, const int iUris); // done + +void uris_meanv(ComMod& com_mod, CmMod& cm_mod, const int iUris); // done + +void uris_update_disp(ComMod& com_mod, CmMod& cm_mod); + +void uris_find_tetra(ComMod& com_mod, CmMod& cm_mod, const int iUris); + +void inside_tet(ComMod& com_mod, int& eNoN, Vector& xp, + Array& xl, int& flag, bool ext); // done + +void uris_read_msh(Simulation* simulation); // done + +void uris_write_vtus(ComMod& com_mod); // done + +void uris_calc_sdf(ComMod& com_mod); // done + +void uris_read_sv(Simulation* simulation, mshType& mesh, const URISFaceParameters* mesh_param); //done + +int in_poly(Vector& P, Array& P1, bool ext); // done + +int same_side(Vector& v1, Vector& v2, Vector& v3, + Vector& v4, Vector& p, bool ext); //done + +} + +#endif + diff --git a/Code/Source/solver/utils.cpp b/Code/Source/solver/utils.cpp index 690aa15f6..cb9d9633c 100644 --- a/Code/Source/solver/utils.cpp +++ b/Code/Source/solver/utils.cpp @@ -399,4 +399,18 @@ void swap(int& value1, int& value2) std::swap(value1, value2); } -}; +/// @brief Reproduces 'FINDLOC' in Fortran for 2D arrays. +void find_loc(const Array& array, int value, std::array& ind) +{ + ind = {-1, -1}; + for (size_t i = 0; i < array.nrows(); ++i) { + for (size_t j = 0; j < array.ncols(); ++j) { + if (array(i,j) == value) { + ind = {static_cast(i), static_cast(j)}; + return; // Exit the function after finding the first occurrence + } + } + } +} + +}; \ No newline at end of file diff --git a/Code/Source/solver/utils.h b/Code/Source/solver/utils.h index 515c4a4bf..867d90dd9 100644 --- a/Code/Source/solver/utils.h +++ b/Code/Source/solver/utils.h @@ -89,6 +89,8 @@ int sign(double value); void swap(int& value1, int& value2); +void find_loc(const Array& array, int value, std::array& ind); + }; #endif diff --git a/Code/Source/solver/vtk_xml.cpp b/Code/Source/solver/vtk_xml.cpp index 256e31318..80721a46a 100644 --- a/Code/Source/solver/vtk_xml.cpp +++ b/Code/Source/solver/vtk_xml.cpp @@ -1000,6 +1000,12 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array outNames(nOut); std::vector outS(nOut+1); std::vectoroutNamesE(nOute); @@ -1312,6 +1318,22 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array(com_mod.uris[iUris].sdf(Ac)); + } + } + } + } // iM for loop @@ -1392,6 +1414,8 @@ void write_vtus(Simulation* simulation, const Array& lA, const Array& lA, const Arrayset_point_data(outNames[iOut], tmpV); } @@ -1437,7 +1464,6 @@ void write_vtus(Simulation* simulation, const Array& lA, const Arrayset_element_data("Domain_ID", tmpI); } - if (!com_mod.savedOnce) { com_mod.savedOnce = true; ne = ne + 1; @@ -1454,7 +1480,6 @@ void write_vtus(Simulation* simulation, const Array& lA, const Arrayset_element_data("Proc_ID", tmpI); } } - // Write the mesh ID // if (nMsh > 1) { @@ -1468,7 +1493,6 @@ void write_vtus(Simulation* simulation, const Array& lA, const Arrayset_element_data("Mesh_ID", tmpI); } } // if (com_mod.savedOnce || nMsh > 1) - // Write element Jacobian and von Mises stress if necessary // for (int l = 0; l < nOute; l++) { @@ -1486,7 +1510,6 @@ void write_vtus(Simulation* simulation, const Array& lA, const Arrayset_element_data(outNamesE[l], tmpVe); } - // Write element ghost cells if necessary if (lIbl) { ne = ne + 1; @@ -1502,7 +1525,6 @@ void write_vtus(Simulation* simulation, const Array& lA, const Arrayset_element_data("EGHOST", tmpI); } - vtk_writer->write(); delete vtk_writer; } diff --git a/tests/cases/ris/pipe_ris_3d/README.md b/tests/cases/ris/pipe_ris_3d/README.md new file mode 100644 index 000000000..b89bb1b93 --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/README.md @@ -0,0 +1,16 @@ +# **Problem Description** +Simulation of a 3D pipe flow problem with an explicit RIS (resistive immersed surfaces) valve model. + + + +The model parameters are specified in the `Add_RIS_projection` sub-section +``` + + right_ris + 1.e6 + 1.e-8 + +``` + +## Reference +Astorino, M., Hamers, J., Shadden, S. C., & Gerbeau, J. F. (2012). A robust and efficient valve model based on resistive immersed surfaces. International Journal for Numerical Methods in Biomedical Engineering, 28(9), 937-959. https://doi.org/10.1002/cnm.2474. \ No newline at end of file diff --git a/tests/cases/ris/pipe_ris_3d/inlet.dat b/tests/cases/ris/pipe_ris_3d/inlet.dat new file mode 100644 index 000000000..197a968eb --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/inlet.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:55fd17868e0fbc70b44f079bbffb45f9d47d022b754a4a2c5ebf4f3c7c20f85f +size 26507 diff --git a/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-complete.mesh.vtu b/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-complete.mesh.vtu new file mode 100644 index 000000000..b1695d750 --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7b706a0115ee1f17167b4ffdcb2673c1d554685b2fe7cda5823ec46219b17943 +size 82271 diff --git a/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_1.vtp b/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_1.vtp new file mode 100644 index 000000000..56144cdcd --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5c6bbf4ecbd6f3700d817baaa68bf2ba3ff9b24eb7d8dd8059487d12221cdc91 +size 29841 diff --git a/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_2.vtp b/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_2.vtp new file mode 100644 index 000000000..e619aee2e --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_2.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cafd1635d10ce0fd061c89924c350651d34000c2ce74e7a68ee97399bb3d9648 +size 5850 diff --git a/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_3.vtp b/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_3.vtp new file mode 100644 index 000000000..40e95d560 --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/mesh/1-mesh-complete/mesh-surfaces/noname_3.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8d78ecfbb274e0cd125179de062122373fa3fa95f7f33dd5809b6e48748c558f +size 5816 diff --git a/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-complete.mesh.vtu b/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-complete.mesh.vtu new file mode 100644 index 000000000..b09bc9afe --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:79bd6fbe4f3dfed7e36a6761da27112f76ba6bedefda051c8c92b607345ea58d +size 82585 diff --git a/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_1.vtp b/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_1.vtp new file mode 100644 index 000000000..94d816087 --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:241f93464507eccaaf806c5963f2b90ddda432093264dd8a7fa645e183d5ecd9 +size 5840 diff --git a/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_2.vtp b/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_2.vtp new file mode 100644 index 000000000..b4898a731 --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_2.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9bba18a86cea41bc677d1aae1b5931116cc99835035ddc2b7af1e37fffb2d0bd +size 29794 diff --git a/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_3.vtp b/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_3.vtp new file mode 100644 index 000000000..ed5431342 --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/mesh/2-mesh-complete/mesh-surfaces/noname_3.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:54399c0a653a9866a7b570a3d12216789acdf75a6ac9ad6b87b61fd52434de7b +size 5858 diff --git a/tests/cases/ris/pipe_ris_3d/result_005.vtu b/tests/cases/ris/pipe_ris_3d/result_005.vtu new file mode 100644 index 000000000..bff53985b --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/result_005.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b3a0dc827241c862fcc23d5ab30f6039e76300751f64458ad02d218a1ec1a279 +size 439958 diff --git a/tests/cases/ris/pipe_ris_3d/solver.xml b/tests/cases/ris/pipe_ris_3d/solver.xml new file mode 100644 index 000000000..1e51a0e90 --- /dev/null +++ b/tests/cases/ris/pipe_ris_3d/solver.xml @@ -0,0 +1,200 @@ + + + + + + true + 3 + 5 + 0.001 + 0.50 + STOP_SIM + + 1 + result + 5 + 1 + + 10 + 0 + + 1 + 0 + 1 + + + + + + mesh/1-mesh-complete/mesh-complete.mesh.vtu + + mesh/1-mesh-complete/mesh-surfaces/noname_3.vtp + + + mesh/1-mesh-complete/mesh-surfaces/noname_2.vtp + + + mesh/1-mesh-complete/mesh-surfaces/noname_1.vtp + + 0 + + + + mesh/2-mesh-complete/mesh-complete.mesh.vtu + + mesh/2-mesh-complete/mesh-surfaces/noname_1.vtp + + + mesh/2-mesh-complete/mesh-surfaces/noname_3.vtp + + + mesh/2-mesh-complete/mesh-surfaces/noname_2.vtp + + 1 + + + + + + right_ris + 1.e6 + 1.e-8 + + + + + true + 3 + 10 + 1e-3 + + + fluid + 1.06 + + 0.035 + + 1. + + + + fluid + 1.06 + + 0.035 + + 1. + + + + true + true + true + true + + + + true + + + + true + true + + + + + fsils + + 1e-3 + 30 + 50 + + + + Neu + Unsteady + inlet.dat + Flat + true + + + + Neu + Steady + true + 0.0 + + + + Dir + Steady + 0.0 + + + + Dir + Steady + 0.0 + + + + + + + true + 1 + 5 + 1e-6 + 0. + + + + fsils + + 1e-4 + + + + true + + + + Dir + Steady + 0.0 + + + + Dir + Steady + 0.0 + + + + Dir + Steady + 0.0 + + + + Dir + Steady + 0.0 + + + + Dir + Steady + 0.0 + + + + Dir + Steady + 0.0 + + + + + \ No newline at end of file diff --git a/tests/cases/uris/pipe_uris_3d/README.md b/tests/cases/uris/pipe_uris_3d/README.md new file mode 100644 index 000000000..8ce96abf6 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/README.md @@ -0,0 +1,18 @@ +# **Problem Description** +Simulation of a 3D pipe flow problem with an unfitted RIS (resistive immersed surfaces) valve model. + + + +The model parameters are specified in the `Add_URIS_mesh` sub-section +``` + + + meshes/LCC_mesh.vtu + meshes/LCC_motion_open.dat + meshes/LCC_motion_close.dat + + 1.0 + 0.25 + meshes/normal.dat + +``` \ No newline at end of file diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25.vtu b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25.vtu new file mode 100644 index 000000000..a14631021 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7feba70a7c351d8c48d8a5d778d19cd7f6e5339f6f7f5b7d632c1f5f2ebe2094 +size 14475 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25motion_close.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25motion_close.dat new file mode 100644 index 000000000..615f51dea --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25motion_close.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:581461715029ba614dbc802150eb9fb83df26268d824383b840c8e0330235701 +size 226254 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25motion_open.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25motion_open.dat new file mode 100644 index 000000000..5243ff670 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed+25motion_open.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6cdad75be7fe50dba0254a2fa62d3fb94e59524abd28bdce911ef51539ae42c6 +size 226279 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25.vtu b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25.vtu new file mode 100644 index 000000000..e0a5ddd9f --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7d8dbac27e2ce90397508a566910c36208688790216b53ce6018eab527e2c7fa +size 14498 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25motion_close.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25motion_close.dat new file mode 100644 index 000000000..66f38a9ca --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25motion_close.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:30427f3c6ac9f2084060b3d910f7c16b238d9dc8870707378c2a3295a15c121f +size 229256 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25motion_open.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25motion_open.dat new file mode 100644 index 000000000..a59354f4f --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_LCC_morphed-25motion_open.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5e1d4f5d409353cc8066491fdb6904168a8edbc0ab30b808ef06af47d438b437 +size 229282 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25.vtu b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25.vtu new file mode 100644 index 000000000..377bca42d --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2b4091d7f776de027e84e38dbf6cefe3856c096e193da426b9e3559c9d5c51bc +size 13111 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25motion_close.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25motion_close.dat new file mode 100644 index 000000000..3bb69192e --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25motion_close.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ce78f987dbd863dc2832e6ee0d807aa7ed1ccf6767c52b78df366dbde590acfc +size 202352 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25motion_open.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25motion_open.dat new file mode 100644 index 000000000..b4f31164b --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed+25motion_open.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:819ca1f18177653708c8adbd5b921f7f763a9acc46713769c45613f413c34a3e +size 202359 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25.vtu b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25.vtu new file mode 100644 index 000000000..19ad91d52 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5871abaac13b60eb668be75e279b3883714c7c66381c0c6e17bf63be7d251c9f +size 13107 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25motion_close.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25motion_close.dat new file mode 100644 index 000000000..7e5ff6263 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25motion_close.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:133044246c7444ae983f112e6505ddc8dca87f8c6210a6bae625ba8538d88ca2 +size 205014 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25motion_open.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25motion_open.dat new file mode 100644 index 000000000..4f1d89b0b --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_NCC_morphed-25motion_open.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f733e3d7185240f3d361c4a7cde034e183b4002588d6568d1a62d37ffe431ade +size 205021 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25.vtu b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25.vtu new file mode 100644 index 000000000..5d04c8d2f --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:322d7b952218f6b42428b7c0519d0304ddeb89086b5d942115a8005d2fe7af06 +size 14749 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25motion_close.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25motion_close.dat new file mode 100644 index 000000000..41a57fbab --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25motion_close.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9bd17cd1fd169553aca9c923f778c8fc852afa32b3f6586c26807030ed565b66 +size 234154 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25motion_open.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25motion_open.dat new file mode 100644 index 000000000..5e1d76dad --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed+25motion_open.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c56893d8c482bb2d4dd54df8d4cfb1c6cff9b3171038bab3d25ca6b4889149f8 +size 234145 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25.vtu b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25.vtu new file mode 100644 index 000000000..4aaa7444c --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:80f8d9d5f5c0d27ba894800697d982b1b1bcde2fd31775b379f8c135dab462af +size 14729 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25motion_close.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25motion_close.dat new file mode 100644 index 000000000..34ae92af4 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25motion_close.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fbaf9e52810299e9b2839d48de1b07d726ce1dbbdea66c2325c2f9ad8406a1a1 +size 237212 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25motion_open.dat b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25motion_open.dat new file mode 100644 index 000000000..3cbf647b1 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/PAT003_07_RCC_morphed-25motion_open.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8c079625d908daa595e827baf75b79ddec972dd2a085381c950c0533894d92e2 +size 237203 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-complete.mesh.vtu b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-complete.mesh.vtu new file mode 100644 index 000000000..50dadfe2c --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0f65577877cc38f3538e2cf25c57c0ae4ad0265b46bbf531acbe1f97eb6a0fe2 +size 230860 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_1.vtp b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_1.vtp new file mode 100644 index 000000000..0efdaff6c --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:aded76c4d51fe9ae782ff64e288602f16eee6bb6333848af783d93f4635ffff0 +size 37570 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_2.vtp b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_2.vtp new file mode 100644 index 000000000..0dcf01c24 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_2.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0e2c6096bf67fb875c6ca6c72b6e01e01056b8965552b4f3bb39e24fa10d5abf +size 5302 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_3.vtp b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_3.vtp new file mode 100644 index 000000000..83d0b65ef --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/mesh-surfaces/noname_3.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f35edcf3d495bd884df8eca59891d2baba87ca39217bae9e78237338da5062d5 +size 5397 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/1_displacement.dat b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/1_displacement.dat new file mode 100644 index 000000000..0468052c1 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/1_displacement.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:74cd9de4a9d704bc93f4b498c3b4c658328980758d3b2e3f52f04a2d00262323 +size 51549656 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/2_displacement.dat b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/2_displacement.dat new file mode 100644 index 000000000..cb8e2e5ea --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/2_displacement.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7cceb6ce23e9415af96da385c489894f6a4a970accf310f8878118def57dd462 +size 4262556 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/3_displacement.dat b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/3_displacement.dat new file mode 100644 index 000000000..ee13562e2 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/3_displacement.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:57c2021a5dfa07111b82c9da56ff91dde6651c9fadf29e8574fc0f408f48297f +size 4475766 diff --git a/tests/cases/uris/pipe_uris_3d/meshes/normal.dat b/tests/cases/uris/pipe_uris_3d/meshes/normal.dat new file mode 100644 index 000000000..0a0f7d584 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/meshes/normal.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0eeca5df6aff4a09db5908e091542decb1228d17f6c42425bcf0c48bf169c751 +size 48 diff --git a/tests/cases/uris/pipe_uris_3d/result_005.vtu b/tests/cases/uris/pipe_uris_3d/result_005.vtu new file mode 100644 index 000000000..025a881eb --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/result_005.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:884dedfcfb99696dfac3f47521f52b3eea2663d9b22a7564eddce4a9035ab44b +size 734173 diff --git a/tests/cases/uris/pipe_uris_3d/solver.xml b/tests/cases/uris/pipe_uris_3d/solver.xml new file mode 100644 index 000000000..ccd61e375 --- /dev/null +++ b/tests/cases/uris/pipe_uris_3d/solver.xml @@ -0,0 +1,196 @@ + + + + + + false + 3 + 5 + 0.001 + 0.20 + STOP_SIM + + 1 + result + 1 + 1 + 5 + 0 + false + 0 + + 1 + 1 + 1 + + + + + + + + meshes/cylinder-mesh-complete/mesh-complete.mesh.vtu + + meshes/cylinder-mesh-complete/mesh-surfaces/noname_1.vtp + + + meshes/cylinder-mesh-complete/mesh-surfaces/noname_3.vtp + + + meshes/cylinder-mesh-complete/mesh-surfaces/noname_2.vtp + + 0 + + + + + + meshes/PAT003_07_LCC_morphed-25.vtu + meshes/PAT003_07_LCC_morphed-25motion_open.dat + meshes/PAT003_07_LCC_morphed-25motion_close.dat + + + meshes/PAT003_07_NCC_morphed-25.vtu + meshes/PAT003_07_NCC_morphed-25motion_open.dat + meshes/PAT003_07_NCC_morphed-25motion_close.dat + + + meshes/PAT003_07_RCC_morphed-25.vtu + meshes/PAT003_07_RCC_morphed-25motion_open.dat + meshes/PAT003_07_RCC_morphed-25motion_close.dat + + 1.0 + 0.25 + meshes/normal.dat + + + + + meshes/PAT003_07_LCC_morphed+25.vtu + meshes/PAT003_07_LCC_morphed+25motion_open.dat + meshes/PAT003_07_LCC_morphed+25motion_close.dat + + + meshes/PAT003_07_NCC_morphed+25.vtu + meshes/PAT003_07_NCC_morphed+25motion_open.dat + meshes/PAT003_07_NCC_morphed+25motion_close.dat + + + meshes/PAT003_07_RCC_morphed+25.vtu + meshes/PAT003_07_RCC_morphed+25motion_open.dat + meshes/PAT003_07_RCC_morphed+25motion_close.dat + + 1.0 + 0.25 + meshes/normal.dat + + + + + + true + 1 + 10 + 1e-4 + + + fluid + 1.06 + + 0.04 + + 0.3 + + + + + fsils + + 1e-4 + 1e-4 + 100 + 50 + + + + false + true + true + true + true + true + + + + Neu + Steady + + + + Neu + Steady + + + + Dir + General + meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/1_displacement.dat + Flat + 0 + 0 + 1 + + + + + + true + 1 + 10 + 1e-6 + 0. + + + + fsils + + 1e-4 + + + + true + + + + Dir + General + meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/1_displacement.dat + Flat + 0 + 0 + 1 + + + + Dir + General + meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/3_displacement.dat + Flat + 0 + 0 + 1 + + + + Dir + General + meshes/cylinder-mesh-complete/motion/motion_awave/mesh-complete/2_displacement.dat + Flat + 0 + 0 + 1 + + + + + \ No newline at end of file diff --git a/tests/test_ris.py b/tests/test_ris.py new file mode 100644 index 000000000..b4d3ce90a --- /dev/null +++ b/tests/test_ris.py @@ -0,0 +1,14 @@ +from .conftest import run_with_reference +import os +import subprocess + +# Common folder for all tests in this file +base_folder = "ris" + +# Fields to test +fields = ["Displacement", "Pressure", "Velocity", "Traction", "WSS"] + +def test_pipe_ris_3d(n_proc): + test_folder = "pipe_ris_3d" + t_max = 5 + run_with_reference(base_folder, test_folder, fields, n_proc, t_max) \ No newline at end of file diff --git a/tests/test_uris.py b/tests/test_uris.py new file mode 100644 index 000000000..f33f1817a --- /dev/null +++ b/tests/test_uris.py @@ -0,0 +1,14 @@ +from .conftest import run_with_reference +import os +import subprocess + +# Common folder for all tests in this file +base_folder = "uris" + +# Fields to test +fields = ["Displacement", "Pressure", "Velocity", "WSS"] + +def test_pipe_uris_3d(n_proc): + test_folder = "pipe_uris_3d" + t_max = 5 + run_with_reference(base_folder, test_folder, fields, n_proc, t_max) \ No newline at end of file