/** * \file * \copyright * Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) * Distributed under a Modified BSD License. * See accompanying file LICENSE.txt or * http://www.opengeosys.org/project/license * */ #include "ProcessVariable.h" #include #include #include "BaseLib/Logging.h" #include "BaseLib/Algorithm.h" #include "BaseLib/TimeInterval.h" #include "MeshGeoToolsLib/ConstructMeshesFromGeometries.h" #include "MeshLib/Mesh.h" #include "MeshLib/Node.h" #include "ParameterLib/Utils.h" #include "ProcessLib/BoundaryCondition/BoundaryCondition.h" #include "ProcessLib/BoundaryCondition/CreateBoundaryCondition.h" #include "ProcessLib/BoundaryCondition/DirichletBoundaryConditionWithinTimeInterval.h" #include "ProcessLib/SourceTerms/CreateSourceTerm.h" #include "ProcessLib/SourceTerms/SourceTerm.h" namespace { MeshLib::Mesh const& findMeshInConfig( BaseLib::ConfigTree const& config, std::vector> const& meshes) { // // Get the mesh name from the config. // std::string mesh_name; // Either given directly in or constructed // from _. #ifdef DOXYGEN_DOCU_ONLY //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__mesh} config.getConfigParameterOptional("mesh"); #endif // DOXYGEN_DOCU_ONLY auto optional_mesh_name = //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__mesh} config.getConfigParameterOptional("mesh"); if (optional_mesh_name) { mesh_name = *optional_mesh_name; } else { #ifdef DOXYGEN_DOCU_ONLY //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__geometrical_set} config.getConfigParameterOptional("geometrical_set"); //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__geometry} config.getConfigParameter("geometry"); #endif // DOXYGEN_DOCU_ONLY // Looking for the mesh created before for the given geometry. auto const geometrical_set_name = //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__geometrical_set} config.getConfigParameter("geometrical_set"); auto const geometry_name = //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__geometry} config.getConfigParameter("geometry"); mesh_name = MeshGeoToolsLib::meshNameFromGeometry(geometrical_set_name, geometry_name); } // // Find and extract mesh from the list of meshes. // auto const& mesh = *BaseLib::findElementOrError( begin(meshes), end(meshes), [&mesh_name](auto const& mesh) { assert(mesh != nullptr); return mesh->getName() == mesh_name; }, "Required mesh with name '" + mesh_name + "' not found."); DBUG("Found mesh '{:s}' with id {:d}.", mesh.getName(), mesh.getID()); return mesh; } } // namespace namespace ProcessLib { ProcessVariable::ProcessVariable( BaseLib::ConfigTree const& config, MeshLib::Mesh& mesh, std::vector> const& meshes, std::vector> const& parameters) : //! \ogs_file_param{prj__process_variables__process_variable__name} _name(config.getConfigParameter("name")), _mesh(mesh), //! \ogs_file_param{prj__process_variables__process_variable__components} _n_components(config.getConfigParameter("components")), //! \ogs_file_param{prj__process_variables__process_variable__order} _shapefunction_order(config.getConfigParameter("order")), _deactivated_subdomains(createDeactivatedSubdomains(config, mesh)), _initial_condition(ParameterLib::findParameter( //! \ogs_file_param{prj__process_variables__process_variable__initial_condition} config.getConfigParameter("initial_condition"), parameters, _n_components, &mesh)) { DBUG("Constructing process variable {:s}", _name); if (_shapefunction_order < 1 || 2 < _shapefunction_order) { OGS_FATAL("The given shape function order {:d} is not supported", _shapefunction_order); } // Boundary conditions if (auto bcs_config = //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions} config.getConfigSubtreeOptional("boundary_conditions")) { for ( auto bc_config : //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition} bcs_config->getConfigSubtreeList("boundary_condition")) { auto const& mesh = findMeshInConfig(bc_config, meshes); auto component_id = //! \ogs_file_param{prj__process_variables__process_variable__boundary_conditions__boundary_condition__component} bc_config.getConfigParameterOptional("component"); if (!component_id && _n_components == 1) { // default value for single component vars. component_id = 0; } _bc_configs.emplace_back(std::move(bc_config), mesh, component_id); } } else { INFO("No boundary conditions for process variable '{:s}' found.", _name); } // Source terms //! \ogs_file_param{prj__process_variables__process_variable__source_terms} if (auto sts_config = config.getConfigSubtreeOptional("source_terms")) { for ( auto st_config : //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term} sts_config->getConfigSubtreeList("source_term")) { MeshLib::Mesh const& mesh = findMeshInConfig(st_config, meshes); auto component_id = //! \ogs_file_param{prj__process_variables__process_variable__source_terms__source_term__component} st_config.getConfigParameterOptional("component"); if (!component_id && _n_components == 1) { // default value for single component vars. component_id = 0; } _source_term_configs.emplace_back(std::move(st_config), mesh, component_id); } } else { INFO("No source terms for process variable '{:s}' found.", _name); } } ProcessVariable::ProcessVariable(ProcessVariable&& other) : _name(std::move(other._name)), _mesh(other._mesh), _n_components(other._n_components), _shapefunction_order(other._shapefunction_order), _deactivated_subdomains(std::move(other._deactivated_subdomains)), _initial_condition(std::move(other._initial_condition)), _bc_configs(std::move(other._bc_configs)), _source_term_configs(std::move(other._source_term_configs)) { } std::string const& ProcessVariable::getName() const { return _name; } MeshLib::Mesh const& ProcessVariable::getMesh() const { return _mesh; } std::vector> ProcessVariable::createBoundaryConditions( const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id, unsigned const integration_order, std::vector> const& parameters, Process const& process) { std::vector> bcs; bcs.reserve(_bc_configs.size()); for (auto& config : _bc_configs) { auto bc = createBoundaryCondition( config, dof_table, _mesh, variable_id, integration_order, _shapefunction_order, parameters, process); #ifdef USE_PETSC if (bc == nullptr) { continue; } #endif // USE_PETSC bcs.push_back(std::move(bc)); } if (_deactivated_subdomains.empty()) { return bcs; } createBoundaryConditionsForDeactivatedSubDomains(dof_table, variable_id, parameters, bcs); return bcs; } void ProcessVariable::createBoundaryConditionsForDeactivatedSubDomains( const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id, std::vector> const& parameters, std::vector>& bcs) { auto& parameter = ParameterLib::findParameter( DeactivatedSubdomain::zero_parameter_name, parameters, 1); for (auto const& deactivated_subdomain : _deactivated_subdomains) { auto const& deactivated_subdomain_meshes = deactivated_subdomain->deactivated_subdomain_meshes; for (auto const& deactivated_subdomain_mesh : deactivated_subdomain_meshes) { for (int component_id = 0; component_id < dof_table.getNumberOfVariableComponents(variable_id); component_id++) { // Copy the time interval. std::unique_ptr time_interval = std::make_unique( *deactivated_subdomain->time_interval); auto bc = std::make_unique< DirichletBoundaryConditionWithinTimeInterval>( std::move(time_interval), parameter, *(deactivated_subdomain_mesh->mesh), deactivated_subdomain_mesh->inactive_nodes, dof_table, variable_id, component_id); #ifdef USE_PETSC // TODO: make it work under PETSc too. if (bc == nullptr) { continue; } #endif // USE_PETSC bcs.push_back(std::move(bc)); } } } } void ProcessVariable::updateDeactivatedSubdomains(double const time) { if (_deactivated_subdomains.empty()) { _ids_of_active_elements.clear(); return; } auto found_a_set = std::find_if(_deactivated_subdomains.begin(), _deactivated_subdomains.end(), [&](auto& _deactivated_subdomain) { return _deactivated_subdomain->includesTimeOf(time); }); if (found_a_set == _deactivated_subdomains.end()) { _ids_of_active_elements.clear(); return; } // Already initialized. if (!_ids_of_active_elements.empty()) { return; } auto const& deactivated_materialIDs = (*found_a_set)->materialIDs; auto const* const material_ids = MeshLib::materialIDs(_mesh); _ids_of_active_elements.clear(); auto const number_of_elements = _mesh.getNumberOfElements(); for (std::size_t i = 0; i < number_of_elements; i++) { if (std::binary_search(deactivated_materialIDs.begin(), deactivated_materialIDs.end(), (*material_ids)[i])) { continue; } _ids_of_active_elements.push_back(_mesh.getElement(i)->getID()); } } std::vector> ProcessVariable::createSourceTerms( const NumLib::LocalToGlobalIndexMap& dof_table, const int variable_id, unsigned const integration_order, std::vector> const& parameters) { std::vector> source_terms; transform(cbegin(_source_term_configs), cend(_source_term_configs), back_inserter(source_terms), [&](auto const& config) { return createSourceTerm(config, dof_table, config.mesh, variable_id, integration_order, _shapefunction_order, parameters); }); return source_terms; } } // namespace ProcessLib