1#ifndef DUNE_COPASI_GRID_MAKE_MULTI_DOMAIN_GRID_HH
2#define DUNE_COPASI_GRID_MAKE_MULTI_DOMAIN_GRID_HH
12#include <dune/grid/common/exceptions.hh>
13#include <dune/grid/common/gridinfo.hh>
14#include <dune/grid/common/mcmgmapper.hh>
15#include <dune/grid/common/rangegenerators.hh>
16#include <dune/grid/io/file/gmshreader.hh>
17#include <dune/grid/uggrid/uggridfactory.hh>
18#include <dune/grid/utility/structuredgridfactory.hh>
20#include <dune/common/parametertree.hh>
24#include <fmt/ranges.h>
43template<Concept::MultiDomainGr
id MDGr
id,
class CellDataType =
double>
44std::tuple<std::unique_ptr<MDGrid>, std::unique_ptr<CellData<typename MDGrid::LevelGridView, CellDataType>>>
47 using HostGrid =
typename MDGrid::HostGrid;
48 using MDEntity =
typename MDGrid::template Codim<0>::Entity;
50 const auto& grid_config = config.sub(
"grid");
51 std::size_t
const dim = grid_config.get(
"dimension", std::size_t{ MDGrid::dimensionworld });
52 if (dim != MDGrid::dimensionworld) {
54 "Executable grid dimension does not match with input grid dimension!");
58 auto& compartments_config = config.sub(
"compartments");
60 std::unique_ptr<HostGrid> host_grid_ptr;
61 std::unique_ptr<MDGrid> md_grid_ptr;
62 std::vector<int> gmsh_physical_entity;
64 std::vector<std::pair<std::string,std::function<bool(
const MDEntity&)>>> compartments;
66 if (grid_config.hasKey(
"path")) {
67 auto grid_path = grid_config.template get<std::string>(
"path");
68 spdlog::info(
"Reading grid file '{}'", grid_path);
69 auto host_grid_factory = Dune::GridFactory<HostGrid>{};
70 auto host_grid_parser = Dune::GmshReaderParser<HostGrid>{ host_grid_factory,
false,
false };
72 host_grid_parser.read(grid_path);
73 gmsh_physical_entity = host_grid_parser.elementIndexMap();
75 host_grid_ptr = host_grid_factory.createGrid();
77 std::array<unsigned int, MDGrid::dimensionworld> cells{};
78 std::fill_n(begin(cells), MDGrid::dimensionworld, 1);
79 cells = grid_config.get(
"cells", cells);
81 grid_config.get(
"extensions", FieldVector<double, MDGrid::dimensionworld>(1.));
82 auto lower_left = grid_config.get(
"origin", FieldVector<double, MDGrid::dimensionworld>(0.));
83 upper_right += lower_left;
84 if (MDGrid::dimensionworld == 1) {
86 StructuredGridFactory<HostGrid>::createCubeGrid(lower_left, upper_right, cells);
89 StructuredGridFactory<HostGrid>::createSimplexGrid(lower_left, upper_right, cells);
93 typename MDGrid::MDGridTraits
const md_grid_traits{};
94 md_grid_ptr = std::make_unique<MDGrid>(std::move(host_grid_ptr), md_grid_traits);
96 auto level = grid_config.get(
"refinement_level",
int{ 0 });
98 spdlog::info(
"Applying refinement of level: {}", level);
99 md_grid_ptr->globalRefine(level);
102 std::cout << fmt::format(
"MultiDomainGrid: ");
103 gridinfo(*md_grid_ptr,
" ");
105 auto cell_data = std::make_unique<CellData<typename MDGrid::LevelGridView, CellDataType>>(md_grid_ptr->levelGridView(0));
108 if (not gmsh_physical_entity.empty())
109 cell_data->addData(
"gmsh_id", gmsh_physical_entity);
112 if (grid_config.hasSub(
"cell_data"))
115 std::vector<CellDataType> cell_values(cell_data->size());
116 std::vector<bool> cell_mask(cell_data->size());
119 for (
const auto& compartment : compartments_config.getSubKeys()) {
120 auto& compartment_config = compartments_config.sub(compartment);
121 if (not compartment_config.hasKey(
"id")) {
123 std::size_t
id = compartments.size();
124 compartment_config[
"id"] = std::to_string(
id);
125 const auto& type = compartment_config.get(
"type",
"expression");
126 if (type ==
"expression") {
128 string2parser.at(compartment_config.get(
"parser_type", default_parser_str));
129 std::shared_ptr
const parser_ptr =
make_parser(parser_type);
130 for (std::size_t i = 0; i != cell_data->size(); ++i)
131 parser_ptr->define_variable(cell_data->keys()[i], &cell_values[i]);
132 parser_ptr->set_expression(compartment_config[
"expression"]);
134 parser_context->add_context(*parser_ptr);
135 auto position = std::make_shared<FieldVector<double, MDGrid::dimensionworld>>();
136 for (std::size_t i = 0; i !=
axis_names.size(); ++i) {
137 auto pos_arg = fmt::format(
"position_{}",
axis_names[i]);
138 if (i < MDGrid::dimensionworld) {
139 parser_ptr->define_variable(pos_arg, &(*position)[i]);
141 parser_ptr->define_constant(pos_arg, 0.);
144 parser_ptr->compile();
145 compartments.emplace_back(compartment, [position, parser_ptr, &cell_data, &cell_values, &cell_mask](
const MDEntity& entity) {
146 (*position) = entity.geometry().center();
147 cell_data->getData(entity, cell_values, cell_mask);
148 return FloatCmp::ne(std::invoke(*parser_ptr), 0.);
151 throw format_exception(NotImplemented{},
"Not know type '{}' of compartment", type);
157 if (compartments_config.getSubKeys().size() != compartments.size()) {
158 throw format_exception(InvalidStateException{},
"Compartment ids were setup with wrong number of sub-domains");
160 for (
const auto& compartment : compartments_config.getSubKeys()) {
161 if (not compartments_config.sub(compartment).hasKey(
"id")) {
162 spdlog::warn(
"Compartment '{}' were not assigned an 'id'!", compartment);
166 spdlog::info(
"Applying sub-domain partition");
167 md_grid_ptr->startSubDomainMarking();
168 std::size_t max_domains_per_entity = 0;
169 std::set<std::size_t> used_sub_domains;
170 for (
const auto& cell : elements(md_grid_ptr->leafGridView())) {
171 for (std::size_t
id = max_domains_per_entity = 0;
id != compartments.size(); ++id) {
172 if (compartments[
id].second(cell)) {
173 used_sub_domains.insert(
id);
174 ++max_domains_per_entity;
175 md_grid_ptr->addToSubDomain(
id, cell);
178 if (max_domains_per_entity >
static_cast<std::size_t
>(md_grid_traits.maxSubDomainIndex()) + 1u) {
180 "This version of dune-copasi does not support to"
181 " have more than {} domains per entity!",
182 md_grid_traits.maxSubDomainIndex() + 1);
186 if (used_sub_domains.size() != compartments.size()) {
187 spdlog::warn(
"Grid partition uses sub-domain{} {} but config file has {} compartment{}",
188 used_sub_domains.size() == 1 ?
"" :
"s", used_sub_domains,
189 compartments.size(), compartments.size() == 1 ?
"" :
"s");
192 md_grid_ptr->preUpdateSubDomains();
193 md_grid_ptr->updateSubDomains();
194 md_grid_ptr->postUpdateSubDomains();
196 for (std::size_t
id = 0;
id != compartments.size(); ++id) {
197 std::cout << fmt::format(
" SubDomain {{{}: {}}}",
id, compartments[
id].first);
198 gridinfo(md_grid_ptr->subDomain(
id),
" ");
201 return std::make_tuple(std::move(md_grid_ptr), std::move(cell_data));
210template<Concept::MultiDomainGr
id MDGr
id>
211std::unique_ptr<MDGrid>
213 std::shared_ptr<const ParserContext> parser_context = {})
215 return std::get<0>(make_multi_domain_grid_with_cell_data<MDGrid>(config, parser_context));
Definition: axis_names.hh:7
std::unique_ptr< Parser > make_parser(ParserType parser_type=default_parser)
auto format_exception(Exception &&e, fmt::format_string< Args... > format, Args &&... args)
Definition: exceptions.hh:23
std::unique_ptr< MDGrid > make_multi_domain_grid(Dune::ParameterTree &config, std::shared_ptr< const ParserContext > parser_context={})
Creates a multi domain grid from a config file From given configuration file with grid and compartmen...
Definition: make_multi_domain_grid.hh:212
auto ostream2spdlog()
Standard output redirection to spdlog.
Definition: ostream_to_spdlog.hh:20
std::vector< std::string > axis_names
void cell_data_parser(const ParameterTree &grid_config, CellData< GV, T > &cell_data)
Parse data file with cell data and add its values to cell data.
Definition: cell_data_parser.hh:94
std::tuple< std::unique_ptr< MDGrid >, std::unique_ptr< CellData< typename MDGrid::LevelGridView, CellDataType > > > make_multi_domain_grid_with_cell_data(Dune::ParameterTree &config, std::shared_ptr< const ParserContext > parser_context={})
Creates a multi domain grid from a config file From given configuration file with grid and compartmen...
Definition: make_multi_domain_grid.hh:45