Dune::Copasi 2.1.0
local_equations.hh
Go to the documentation of this file.
1#ifndef DUNE_COPASI_MODEL_LOCAL_EQUATIONS_LOCAL_EQUATIONS_HH
2#define DUNE_COPASI_MODEL_LOCAL_EQUATIONS_LOCAL_EQUATIONS_HH
3
12
13#include <dune/pdelab/concepts/basis.hh>
14#include <dune/pdelab/common/container_entry.hh>
15
16#include <dune/geometry/dimension.hh>
17
18#include <dune/common/fvector.hh>
19#include <dune/common/indices.hh>
20#include <dune/common/overloadset.hh>
21#include <dune/common/parametertree.hh>
22#include <dune/common/tuplevector.hh>
23
24#include <function2/function2.hpp>
25#include <spdlog/spdlog.h>
26
27#include <functional>
28#include <utility>
29#include <variant>
30#include <set>
31
33
34// this class holds a data-structure for each equation that contains functors to be evaluated.
35// Additionally, it contains the values with respect these functors may be evaluated if they are
36// non-linear All the functors are require to be thread-safe!
37template<std::size_t dim>
38class LocalEquations : public LocalDomain<dim>
39{
40 using Scalar = FieldVector<double, 1>;
41 using Vector = FieldVector<double, dim>;
42 using Tensor = FieldMatrix<double, dim, dim>;
43
44 using CompartmentPath = TypeTree::HybridTreePath<index_constant<0>, std::size_t, std::size_t>;
45 using MembranePath = TypeTree::HybridTreePath<index_constant<1>, std::size_t, std::size_t>;
46
47 LocalEquations() = default;
48public:
51
52 virtual ~LocalEquations() override = default;
53
56
57private:
58 static const Concept::CompartmentScalarLocalBasisNode auto& path_to_local_basis_node(
59 CompartmentPath path,
60 const Concept::CompartmentLocalBasisNode auto& lbasis)
61 {
62 assert(path[Indices::_1] == 0);
63 return lbasis.child(path[Indices::_2]);
64 }
65
66 static const Concept::CompartmentScalarLocalBasisNode auto& path_to_local_basis_node(
67 CompartmentPath path,
69 {
70 return lbasis.child(path[Indices::_1]).child(path[Indices::_2]);
71 }
72
73 static const Concept::MembraneSubEntitiesLocalBasisNode auto& path_to_local_basis_node(
74 MembranePath path,
75 const Concept::MembraneLocalBasisNode auto& lbasis)
76 {
77 assert(path[Indices::_1] == 0);
78 return lbasis.child(path[Indices::_2]);
79 }
80
81 static const Concept::MembraneSubEntitiesLocalBasisNode auto& path_to_local_basis_node(
82 MembranePath path,
83 const Concept::MultiCompartmentLocalBasisNode auto& lbasis)
84 {
85 return lbasis.child(path[Indices::_1]).child(path[Indices::_2]);
86 }
87
88 enum class FactoryFalgs
89 {
90 Reaction = 1 << 0,
91 Diffusion = 1 << 1,
92 Velocity = 1 << 2,
93 Outflow = 1 << 3,
94 Storage = 1 << 4
95 };
96
97public:
98 struct CompartmentNode;
99 struct MembraneNode;
100
101 template<class Signature>
102 struct CompartmentPartialDerivative : public fu2::unique_function<Signature>
103 {
104 CompartmentPartialDerivative(fu2::unique_function<Signature>&& callable,
105 const CompartmentNode& _wrt)
106 : fu2::unique_function<Signature>{ std::move(callable) }
107 , wrt{ _wrt }
108 {
109 }
111 };
112
113 template<class Signature>
114 struct MembranePartialDerivative : public fu2::unique_function<Signature>
115 {
116 MembranePartialDerivative(fu2::unique_function<Signature>&& callable,
117 const MembraneNode& _wrt)
118 : fu2::unique_function<Signature>{ std::move(callable) }
119 , wrt{ _wrt }
120 {
121 }
123 };
124
125 template<class Signature>
126 struct CompartmentDifferentiableFunction : public fu2::unique_function<Signature>
127 {
128 std::vector<CompartmentPartialDerivative<Signature>> compartment_jacobian = {};
129 };
130
131 template<class Signature>
132 struct MembraneDifferentiableFunction : public fu2::unique_function<Signature>
133 {
134 std::vector<CompartmentPartialDerivative<Signature>> compartment_jacobian = {};
135 std::vector<MembranePartialDerivative<Signature>> membrane_jacobian = {};
136 };
137
139 CompartmentDifferentiableFunction<Scalar() const noexcept>;
141 MembraneDifferentiableFunction<Scalar() const noexcept>;
142
144 CompartmentDifferentiableFunction<Vector() const noexcept>;
146 MembraneDifferentiableFunction<Vector() const noexcept>;
147
149 : public CompartmentDifferentiableFunction<Vector(Vector) const noexcept>
150 {
151 CompartmentDiffusionApply(fu2::unique_function<Vector(Vector) const noexcept>&& callable,const CompartmentNode& _wrt)
152 : CompartmentDifferentiableFunction<Vector(Vector) const noexcept>{ std::move(callable) }
153 , wrt{ _wrt }
154 {
155 }
157 };
158
160 : public MembraneDifferentiableFunction<Vector(Vector) const noexcept>
161 {
163 fu2::unique_function<Vector(Vector) const noexcept>&& callable,
164 const MembraneNode& _wrt)
165 : MembraneDifferentiableFunction<Vector(Vector) const noexcept>{ std::move(callable) }
166 , wrt{ _wrt }
167 {
168 }
170 };
171
173 {
174 Scalar& value;
175 Vector& gradient;
176
177 CompartmentPath path;
178 std::string name;
179
183
184 std::vector<CompartmentDiffusionApply> cross_diffusion;
185 std::vector<MembraneScalarFunction> outflow;
186
187 CompartmentNode(Scalar& value, Vector& gradient, CompartmentPath path, const std::string& name) : value{value}, gradient{gradient}, path{path}, name{name} {}
188
190 const PDELab::Concept::LocalBasis auto& lbasis) const
191 {
192 return path_to_local_basis_node(path, lbasis.tree());
193 }
194
195 void debug() const
196 {
197 std::cout << fmt::format("Name: {}\n", name);
198 std::cout << "\tPath: " << path << std::endl;
199 std::cout << fmt::format("\tValue: {}\n", value[0]);
200 std::cout << "\tGradient: " << gradient << std::endl;
201 if (reaction) {
202 std::cout << fmt::format("\tReaction: {}\n", reaction()[0]);
203 for (const auto& jac : reaction.compartment_jacobian)
204 std::cout << fmt::format("\t\tJacobian wrt '{}': {}\n", jac.wrt.name, jac()[0]);
205 }
206 if (storage)
207 std::cout << fmt::format("\tStorage: {}\n", storage()[0]);
208 for (const auto& diffusion : cross_diffusion) {
209 std::cout << fmt::format("\tCross Diffusion wrt '{}': ", diffusion.wrt.name)
210 << diffusion(diffusion.wrt.gradient) << std::endl;
211 for (const auto& jac : diffusion.compartment_jacobian)
212 std::cout << fmt::format("\t\tJacobian wrt '{}': ", jac.wrt.name) << jac(jac.wrt.gradient)
213 << std::endl;
214 }
215 }
216 };
217
219 {
220 Scalar& value;
221 Vector& gradient;
222
223 MembranePath path;
224 std::string name;
225 bool is_linear = false;
226
230
231 std::vector<MembraneDiffusionApply> cross_diffusion;
232 std::vector<MembraneScalarFunction> outflow;
233
235 const PDELab::Concept::LocalBasis auto& lbasis) const
236 {
237 return path_to_local_basis_node(path, lbasis.tree());
238 }
239
240 void debug() const
241 {
242 std::cout << fmt::format("Name: {}\n", name);
243 std::cout << "\tPath: " << path << std::endl;
244 std::cout << fmt::format("\tValue: {}\n", value[0]);
245 std::cout << "\tGradient: " << gradient << std::endl;
246 if (reaction) {
247 std::cout << fmt::format("\tReaction: {}\n", reaction()[0]);
248 for (const auto& jac : reaction.compartment_jacobian)
249 std::cout << fmt::format("\t\tJacobian wrt '{}': {}\n", jac.wrt.name, jac()[0]);
250 }
251 if (storage)
252 std::cout << fmt::format("\tStorage: {}\n", storage()[0]);
253 for (const auto& diffusion : cross_diffusion) {
254 std::cout << fmt::format("\tCross Diffusion wrt '{}': ", diffusion.wrt.name)
255 << diffusion(diffusion.wrt.gradient) << std::endl;
256 for (const auto& jac : diffusion.compartment_jacobian)
257 std::cout << fmt::format("\t\tJacobian wrt '{}': ", jac.wrt.name) << jac(jac.wrt.gradient)
258 << std::endl;
259 }
260 }
261 };
262
263 template<PDELab::Concept::LocalBasis LocalBasis, class CellDataGridView = typename LocalBasis::GlobalBasis::EntitySet, class CellDataType = double>
264 [[nodiscard]] static std::unique_ptr<LocalEquations> make(
265 const LocalBasis& lbasis,
266 const ParameterTree& eqs_cfg = {},
267 const ParameterTree& domain_cfg = {},
268 std::shared_ptr<const FunctorFactory<dim>> functor_factory = nullptr,
269 std::shared_ptr<const CellData<CellDataGridView, CellDataType>> grid_cell_data = nullptr,
271 {
272 TRACE_EVENT("dune", "LocalEquations::make");
273 auto local_values = std::unique_ptr<LocalEquations>(new LocalEquations{});
274
275 // calculate how many nodes there are
276 std::size_t compartment_count = 0;
277 std::size_t membrane_count = 0;
278 PDELab::forEachLeafNode(
279 lbasis.tree(),
280 overload([&](const Concept::CompartmentScalarLocalBasisNode auto&) { ++compartment_count; },
281 [&](const Concept::MembraneScalarLocalBasisNode auto&) { ++membrane_count; }));
282 local_values->_values.reserve(compartment_count + membrane_count);
283 local_values->_gradients.reserve(compartment_count + membrane_count);
284
285 local_values->initialize_nodes(
286 lbasis.tree(), [&](auto path) { return lbasis.globalBasis().subSpace(path).name(); });
287
288 // check that names are not repeated
289 std::set<std::string_view> names;
290 PDELab::forEach(local_values->_nodes, [&](auto& compartments_fncs) {
291 for (const auto& compartment_fncs : compartments_fncs)
292 for (const auto& component_fncs : compartment_fncs)
293 if (auto [it, inserted] = names.insert(component_fncs.name); not inserted)
294 throw format_exception(
295 InvalidStateException{}, "\tVariable with name '{}' is repeated", *it);
296 });
297
298 // configure grid context data
299 if (grid_cell_data) {
300 local_values->cell_values.resize( grid_cell_data->size() );
301 local_values->cell_mask.resize( grid_cell_data->size() );
302 local_values->cell_keys = grid_cell_data->keys();
303 }
304
305 if ((opts.any() or domain_cfg.hasSub("deformation")) and not functor_factory)
306 throw format_exception(InvalidStateException{},
307 "Equations cannot be configured without a functor factory");
308
309 if (opts.any())
310 local_values->configure_eqs(eqs_cfg, functor_factory, opts);
311
312 if (domain_cfg.hasSub("deformation"))
313 set_differentiable_function(*local_values,
314 *functor_factory,
315 local_values->domain_deformation,
316 "domain.deformation",
317 domain_cfg.sub("deformation"),
318 VectorTag{});
319 return local_values;
320 }
321
322 template<class CellDataGridView, class CellDataType = double>
323 static std::unique_ptr<LocalEquations> make_stiffness(
324 const PDELab::Concept::LocalBasis auto& lbasis,
325 const ParameterTree& eqs_cfg,
326 const ParameterTree& domain_cfg,
327 std::shared_ptr<const FunctorFactory<dim>> functor_factory,
328 std::shared_ptr<const CellData<CellDataGridView, CellDataType>> grid_cell_data
329 )
330 {
331 return make(lbasis,
332 eqs_cfg,
333 domain_cfg,
334 functor_factory,
335 grid_cell_data,
336 FactoryFalgs::Reaction | FactoryFalgs::Diffusion | FactoryFalgs::Velocity |
337 FactoryFalgs::Outflow);
338 }
339
340 template<class CellDataGridView, class CellDataType = double>
341 static std::unique_ptr<LocalEquations> make_mass(
342 const PDELab::Concept::LocalBasis auto& lbasis,
343 const ParameterTree& eqs_cfg,
344 const ParameterTree& domain_cfg,
345 std::shared_ptr<const FunctorFactory<dim>> functor_factory,
346 std::shared_ptr<const CellData<CellDataGridView, CellDataType>> grid_cell_data
347 )
348 {
349 return make(lbasis,
350 eqs_cfg,
351 domain_cfg,
352 functor_factory,
353 grid_cell_data,
354 FactoryFalgs::Storage);
355 }
356
357 template<class Tree>
360 const auto& get_equation(const Tree& tree) const
361 {
362 return PDELab::containerEntry(_nodes, make_path(tree));
363 }
364
365 template<class Tree>
368 auto& get_value(const Tree& tree)
369 {
370 return PDELab::containerEntry(_nodes, make_path(tree)).value;
371 }
372
373 template<class Tree>
376 auto& get_gradient(const Tree& tree)
377 {
378 return PDELab::containerEntry(_nodes, make_path(tree)).gradient;
379 }
380
381 void forEachValue(Codim<0>,
382 const std::function<void(std::string_view,
383 const FieldVector<double, 1>&,
384 const FieldVector<double, dim>&)>& apply) const override
385 {
386 PDELab::forEach(nodes(), [&](auto& compartments) {
387 for (auto& compartment_fncs : compartments)
388 for (auto& component_fncs : compartment_fncs)
389 apply(component_fncs.name, component_fncs.value, component_fncs.gradient);
390 });
391 }
392
393 void forEachValue(Codim<1>,
394 const std::function<void(std::string_view,
395 const FieldVector<double, 1>&,
396 const FieldVector<double, dim - 1>&)>&) const override
397 {
398 throw format_exception(NotImplemented{}, "...");
399 }
400
401 void clear()
402 {
403 std::fill(begin(_values), end(_values), 0.);
404 std::fill(begin(_gradients), end(_gradients), 0.);
405 }
406
407 void debug() const
408 {
409 PDELab::forEach(_nodes, [&](auto& compartments) {
410 for (auto& compartment_fncs : compartments)
411 for (auto& func : compartment_fncs)
412 func.debug();
413 });
414 }
415
416 const auto& nodes() const { return _nodes; }
417
419
420private:
421 template<class Tree>
424 auto make_path(const Tree& tree) const
425 {
426 // the tree may be one of many kind of trees
427
428 // membrane/compartment is distinsuihed because the first index is a template parameter
429 // and because they have a different concept requirements
430 auto compartment_type = compartment_index(tree);
431
432 // remove compartment/membrane index (always comes from a tuple and is compile time constant)
433 auto mi = [&]() {
434 auto front_index = front(tree.orderingViewPath());
435 // in some cases, the tree has both comartment types in the tree, then the first index of its
436 // path is an integral constant
437 if constexpr (IsCompileTimeConstant<decltype(front_index)>{}) {
438 assert(front_index == compartment_type);
439 return pop_front(tree.orderingViewPath());
440 } else {
441 // otherwise, the tree has only one type of compartments
442 return tree.orderingViewPath();
443 }
444 }();
445 assert(mi.size() <= 2);
446 // all the compartment trees have their field id as their last path in the ordering view path
447 // (different from tree path if membrane)
448 auto field_id = back(mi);
449 // in case of different compartments, the first index (after stripping its compartment type) is
450 // the compartment id
451 auto compartment_id = (mi.size() == 2) ? front(mi) : 0;
452 return TypeTree::treePath(compartment_type, compartment_id, field_id);
453 }
454
455 // (bulk|membrane, compartment, component)
456 TupleVector<std::vector<std::vector<CompartmentNode>>, std::vector<std::vector<MembraneNode>>>
457 _nodes;
458 std::array<std::vector<std::string>, 2> _compartment_names;
459
460 std::vector<Scalar> _values;
461 std::vector<Vector> _gradients;
462
463 auto compartment_index(const Concept::CompartmentLocalBasisTree auto&) const
464 {
465 return Indices::_0;
466 }
467
468 auto compartment_index(const Concept::MembraneLocalBasisTree auto&) const { return Indices::_1; }
469
470 template<class Tree>
471 requires Concept::CompartmentLocalBasisNode<Tree> || Concept::MembraneLocalBasisNode<Tree>
472 void initialize_nodes(const Tree& tree, auto fname, std::size_t c = 0)
473 {
474 auto ct = compartment_index(tree);
475 _nodes[ct].resize(c + 1);
476 _nodes[ct][c].clear();
477 _compartment_names[ct].resize(c + 1);
478 for (std::size_t i = 0; i != tree.degree(); ++i) {
479 _nodes[ct][c].emplace_back(_values.emplace_back(),
480 _gradients.emplace_back(),
481 TypeTree::treePath(ct, c, i),
482 fname(tree.child(i).orderingViewPath()));
483 _compartment_names[ct][c] = fname(pop_back(tree.child(i).orderingViewPath()));
484 }
485 }
486
487 template<class Tree>
488 requires Concept::MultiCompartmentLocalBasisNode<Tree> ||
489 Concept::MultiMembraneLocalBasisNode<Tree>
490 void initialize_nodes(const Tree& tree, auto fname)
491 {
492 for (std::size_t c = 0; c != tree.degree(); ++c)
493 initialize_nodes(tree.child(c), fname, c);
494 }
495
496 void initialize_nodes(const Concept::MultiCompartmentMembraneLocalBasisNode auto& tree,
497 auto fname)
498 {
499 using namespace Indices;
500 initialize_nodes(tree.child(_0), fname);
501 initialize_nodes(tree.child(_1), fname);
502 }
503
504 void initialize_nodes(const Concept::CompartmentMembraneLocalBasisNode auto& tree, auto fname)
505 {
506 using namespace Indices;
507 initialize_nodes(tree.child(_0), fname, 0);
508 initialize_nodes(tree.child(_1), fname, 0);
509 }
510
511 using ScalarTag = index_constant<0>;
512 using VectorTag = index_constant<1>;
513 using TensorApplyTag = index_constant<2>;
514
515 // scalar functor maker overload
516 static fu2::unique_function<Scalar() const noexcept> make_functor(
517 const LocalEquations& local_eqs,
518 const FunctorFactory<dim>& functor_factory,
519 std::string_view prefix,
520 const ParameterTree& config,
521 int codim,
522 ScalarTag)
523 {
524 return functor_factory.make_scalar(prefix, config, local_eqs, codim);
525 }
526
527 // vector functor maker overload
528 static fu2::unique_function<Vector() const noexcept> make_functor(
529 const LocalEquations& local_eqs,
530 const FunctorFactory<dim>& functor_factory,
531 std::string_view prefix,
532 const ParameterTree& config,
533 int codim,
534 VectorTag)
535 {
536 return functor_factory.make_vector(prefix, config, local_eqs, codim);
537 }
538
539 // tensor functor maker overload
540 static fu2::unique_function<Vector(Vector) const noexcept> make_functor(
541 const LocalEquations& local_eqs,
542 const FunctorFactory<dim>& functor_factory,
543 std::string_view prefix,
544 const ParameterTree& config,
545 int codim,
546 TensorApplyTag)
547 {
548 return functor_factory.make_tensor_apply(prefix, config, local_eqs, codim);
549 }
550
551 // compartment differentiable functor maker overload
552 template<class Signature>
553 static auto set_differentiable_function(const LocalEquations& local_eqs,
554 const FunctorFactory<dim>& functor_factory,
555 CompartmentDifferentiableFunction<Signature>& function,
556 std::string_view prefix,
557 const ParameterTree& function_config,
558 auto range_tag)
559 {
560 function = CompartmentDifferentiableFunction<Signature>{ make_functor(
561 local_eqs, functor_factory, prefix, function_config, 0, range_tag) };
562 std::set<std::string_view> debug_set;
563 const auto& jac_config = function_config.sub("jacobian");
564 for (const auto& compartment_fncs_jac : local_eqs._nodes[Indices::_0])
565 for (const auto& component_fncs_jac : compartment_fncs_jac)
566 if (jac_config.hasSub(component_fncs_jac.name))
567 if (auto jac =
568 make_functor(local_eqs,
569 functor_factory,
570 fmt::format("{}.jacobian.{}", prefix, component_fncs_jac.name),
571 jac_config.sub(component_fncs_jac.name),
572 0,
573 range_tag)) {
574 function.compartment_jacobian.emplace_back(std::move(jac), component_fncs_jac);
575 debug_set.insert(component_fncs_jac.name);
576 }
577 if (jac_config.getSubKeys().size() != debug_set.size())
578 spdlog::warn("Some sub-sections in \"jacobian\" section are being ignored");
579 }
580
581 // membrane differentiable functor maker overload
582 template<class Signature>
583 static auto set_differentiable_function(const LocalEquations& local_eqs,
584 const FunctorFactory<dim>& functor_factory,
585 MembraneDifferentiableFunction<Signature>& function,
586 std::string_view prefix,
587 const ParameterTree& function_config,
588 auto range_tag)
589 {
590 function = MembraneDifferentiableFunction<Signature>{ make_functor(
591 local_eqs, functor_factory, prefix, function_config, 1, range_tag) };
592 const auto& jac_config = function_config.sub("jacobian");
593 for (const auto& compartment_fncs_jac :local_eqs. _nodes[Indices::_0])
594 for (const auto& component_fncs_jac : compartment_fncs_jac)
595 if (jac_config.hasSub(component_fncs_jac.name))
596 if (auto jac =
597 make_functor(local_eqs,
598 functor_factory,
599 fmt::format("{}.jacobian.{}", prefix, component_fncs_jac.name),
600 jac_config.sub(component_fncs_jac.name),
601 1,
602 range_tag))
603 function.compartment_jacobian.emplace_back(std::move(jac), component_fncs_jac);
604 for (const auto& membrane_fncs_jac : local_eqs._nodes[Indices::_1])
605 for (const auto& component_fncs_jac : membrane_fncs_jac)
606 if (jac_config.hasSub(component_fncs_jac.name))
607 if (auto jac =
608 make_functor(local_eqs,
609 functor_factory,
610 fmt::format("{}.jacobian.{}", prefix, component_fncs_jac.name),
611 jac_config.sub(component_fncs_jac.name),
612 1,
613 range_tag))
614 function.membrane_jacobian.emplace_back(std::move(jac), component_fncs_jac);
615 }
616
617 void configure_eqs(const ParameterTree& eqs_cfg,
618 std::shared_ptr<const FunctorFactory<dim>> functor_factory,
619 BitFlags<FactoryFalgs> opts)
620 {
621 PDELab::forEach(
622 _nodes, [&]<class Node>(std::vector<std::vector<Node>>& compartments_fncs, auto l) {
623 for (auto& compartment_fncs : compartments_fncs) {
624 for (Node& component_fncs : compartment_fncs) {
625
626 const auto& component_config = eqs_cfg.sub(component_fncs.name, true);
627
628 if (opts.test(FactoryFalgs::Diffusion) and component_config.hasSub("cross_diffusion")) {
629 const auto& cross_diffusion_config = component_config.sub("cross_diffusion");
630 std::set<std::string_view> debug_set;
631 for (const auto& compartment_fncs_cross_diff : compartments_fncs)
632 for (const auto& component_fncs_cross_diff : compartment_fncs_cross_diff)
633 if (cross_diffusion_config.hasSub(component_fncs_cross_diff.name)) {
634 const auto& diffusion_config =
635 cross_diffusion_config.sub(component_fncs_cross_diff.name);
636 debug_set.insert(component_fncs_cross_diff.name);
637 auto& cross_diffusion = component_fncs.cross_diffusion.emplace_back(
638 nullptr, component_fncs_cross_diff);
639 set_differentiable_function(
640 std::as_const(*this),
641 *functor_factory,
642 cross_diffusion,
643 fmt::format("{}.cross_diffusion", component_fncs.name),
644 diffusion_config,
645 TensorApplyTag{});
646 if (not cross_diffusion)
647 component_fncs.cross_diffusion.pop_back();
648 }
649 if (cross_diffusion_config.getSubKeys().size() != debug_set.size())
650 spdlog::warn(
651 "Some sub-sections in \"{}.cross_diffusion\" section are being ignored",
652 component_fncs.name);
653 }
654
655 if (opts.test(FactoryFalgs::Reaction) and component_config.hasSub("reaction"))
656 set_differentiable_function(std::as_const(*this),
657 *functor_factory,
658 component_fncs.reaction,
659 fmt::format("{}.reaction", component_fncs.name),
660 component_config.sub("reaction"),
661 ScalarTag{});
662
663 if (opts.test(FactoryFalgs::Velocity) and component_config.hasSub("velocity"))
664 set_differentiable_function(std::as_const(*this),
665 *functor_factory,
666 component_fncs.velocity,
667 fmt::format("{}.velocity", component_fncs.name),
668 component_config.sub("velocity"),
669 VectorTag{});
670
671 if (opts.test(FactoryFalgs::Storage) and component_config.hasSub("storage"))
672 set_differentiable_function(std::as_const(*this),
673 *functor_factory,
674 component_fncs.storage,
675 fmt::format("{}.storage", component_fncs.name),
676 component_config.sub("storage"),
677 ScalarTag{});
678
679 if (opts.test(FactoryFalgs::Outflow) and component_config.hasSub("outflow")) {
680 if (l == 1)
681 throw format_exception(NotImplemented{},
682 "Outflow for membranes is not implemented");
683 const auto& boundaries_config = component_config.sub("outflow");
684 for (std::size_t i = 0; i != _compartment_names[l].size(); ++i) {
685 if (boundaries_config.hasSub(_compartment_names[l][i])) {
686 component_fncs.outflow.resize(_compartment_names[l].size());
687 set_differentiable_function(
688 std::as_const(*this),
689 *functor_factory,
690 component_fncs.outflow[i],
691 fmt::format("{}.outflow.{}", component_fncs.name, _compartment_names[l][i]),
692 boundaries_config.sub(_compartment_names[l][i]),
693 ScalarTag{});
694 }
695 }
696 }
697 }
698 }
699 });
700 }
701};
702
703} // namespace Dune::Copasi::DiffusionReaction
704
705#endif // DUNE_COPASI_MODEL_LOCAL_EQUATIONS_LOCAL_EQUATIONS_HH
Bit flags for enumerators.
Definition: bit_flags.hh:47
static constexpr BitFlags< Enum > no_flags()
Bitflag with all flags turned off.
Definition: bit_flags.hh:70
constexpr bool any() const
Checks if any flags are set to true.
Definition: bit_flags.hh:207
Container for cell data of a grid view.
Definition: cell_data.hh:25
Definition: local_equations.hh:39
static std::unique_ptr< LocalEquations > make_stiffness(const PDELab::Concept::LocalBasis auto &lbasis, const ParameterTree &eqs_cfg, const ParameterTree &domain_cfg, std::shared_ptr< const FunctorFactory< dim > > functor_factory, std::shared_ptr< const CellData< CellDataGridView, CellDataType > > grid_cell_data)
Definition: local_equations.hh:323
void forEachValue(Codim< 1 >, const std::function< void(std::string_view, const FieldVector< double, 1 > &, const FieldVector< double, dim - 1 > &)> &) const override
Definition: local_equations.hh:393
LocalEquations & operator=(LocalEquations &&)=delete
void clear()
Definition: local_equations.hh:401
auto & get_gradient(const Tree &tree)
Definition: local_equations.hh:376
const auto & nodes() const
Definition: local_equations.hh:416
LocalEquations(const LocalEquations &)=delete
LocalEquations & operator=(const LocalEquations &)=delete
static std::unique_ptr< LocalEquations > make_mass(const PDELab::Concept::LocalBasis auto &lbasis, const ParameterTree &eqs_cfg, const ParameterTree &domain_cfg, std::shared_ptr< const FunctorFactory< dim > > functor_factory, std::shared_ptr< const CellData< CellDataGridView, CellDataType > > grid_cell_data)
Definition: local_equations.hh:341
CompartmentVectorFunction domain_deformation
Definition: local_equations.hh:418
void forEachValue(Codim< 0 >, const std::function< void(std::string_view, const FieldVector< double, 1 > &, const FieldVector< double, dim > &)> &apply) const override
Definition: local_equations.hh:381
static std::unique_ptr< LocalEquations > make(const LocalBasis &lbasis, const ParameterTree &eqs_cfg={}, const ParameterTree &domain_cfg={}, std::shared_ptr< const FunctorFactory< dim > > functor_factory=nullptr, std::shared_ptr< const CellData< CellDataGridView, CellDataType > > grid_cell_data=nullptr, BitFlags< FactoryFalgs > opts=BitFlags< FactoryFalgs >::no_flags())
Definition: local_equations.hh:264
void debug() const
Definition: local_equations.hh:407
auto & get_value(const Tree &tree)
Definition: local_equations.hh:368
const auto & get_equation(const Tree &tree) const
Definition: local_equations.hh:360
Definition: functor_factory.hh:24
Definition: functor_factory_parser.impl.hh:22
Definition: compartment_tree.hh:30
Definition: compartment_tree.hh:86
Definition: compartment_tree.hh:36
Definition: compartment_tree.hh:92
Definition: compartment_tree.hh:25
Definition: factory.hh:28
auto format_exception(Exception &&e, fmt::format_string< Args... > format, Args &&... args)
Definition: exceptions.hh:23
std::vector< CompartmentPartialDerivative< Signature > > compartment_jacobian
Definition: local_equations.hh:128
const CompartmentNode & wrt
Definition: local_equations.hh:156
CompartmentDiffusionApply(fu2::unique_function< Vector(Vector) const noexcept > &&callable, const CompartmentNode &_wrt)
Definition: local_equations.hh:151
CompartmentScalarFunction storage
Definition: local_equations.hh:181
std::vector< MembraneScalarFunction > outflow
Definition: local_equations.hh:185
std::vector< CompartmentDiffusionApply > cross_diffusion
Definition: local_equations.hh:184
std::string name
Definition: local_equations.hh:178
void debug() const
Definition: local_equations.hh:195
CompartmentScalarFunction reaction
Definition: local_equations.hh:180
CompartmentNode(Scalar &value, Vector &gradient, CompartmentPath path, const std::string &name)
Definition: local_equations.hh:187
Vector & gradient
Definition: local_equations.hh:175
const Concept::CompartmentScalarLocalBasisNode auto & to_local_basis_node(const PDELab::Concept::LocalBasis auto &lbasis) const
Definition: local_equations.hh:189
CompartmentPath path
Definition: local_equations.hh:177
Scalar & value
Definition: local_equations.hh:174
CompartmentVectorFunction velocity
Definition: local_equations.hh:182
const CompartmentNode & wrt
Definition: local_equations.hh:110
CompartmentPartialDerivative(fu2::unique_function< Signature > &&callable, const CompartmentNode &_wrt)
Definition: local_equations.hh:104
std::vector< MembranePartialDerivative< Signature > > membrane_jacobian
Definition: local_equations.hh:135
std::vector< CompartmentPartialDerivative< Signature > > compartment_jacobian
Definition: local_equations.hh:134
const MembraneNode & wrt
Definition: local_equations.hh:169
MembraneDiffusionApply(fu2::unique_function< Vector(Vector) const noexcept > &&callable, const MembraneNode &_wrt)
Definition: local_equations.hh:162
void debug() const
Definition: local_equations.hh:240
MembranePath path
Definition: local_equations.hh:223
Scalar & value
Definition: local_equations.hh:220
std::vector< MembraneDiffusionApply > cross_diffusion
Definition: local_equations.hh:231
bool is_linear
Definition: local_equations.hh:225
MembraneScalarFunction storage
Definition: local_equations.hh:228
std::string name
Definition: local_equations.hh:224
Vector & gradient
Definition: local_equations.hh:221
const Concept::MembraneSubEntitiesLocalBasisNode auto & to_local_basis_node(const PDELab::Concept::LocalBasis auto &lbasis) const
Definition: local_equations.hh:234
MembraneScalarFunction reaction
Definition: local_equations.hh:227
std::vector< MembraneScalarFunction > outflow
Definition: local_equations.hh:232
MembraneVectorFunction velocity
Definition: local_equations.hh:229
const MembraneNode & wrt
Definition: local_equations.hh:122
MembranePartialDerivative(fu2::unique_function< Signature > &&callable, const MembraneNode &_wrt)
Definition: local_equations.hh:116
Definition: local_domain.hh:15