Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #ifndef LIBMESH_COMPOSITE_FUNCTION_H
19 : #define LIBMESH_COMPOSITE_FUNCTION_H
20 :
21 : // libMesh includes
22 : #include "libmesh/dense_vector.h"
23 : #include "libmesh/function_base.h"
24 : #include "libmesh/int_range.h"
25 : #include "libmesh/libmesh.h"
26 : #include "libmesh/point.h"
27 :
28 : // C++ includes
29 : #include <algorithm>
30 : #include <utility>
31 : #include <vector>
32 :
33 : namespace libMesh
34 : {
35 :
36 : /**
37 : * \brief A function that returns a vector whose components are defined by
38 : * multiple functions.
39 : *
40 : * A function which is defined by composing the result of different functions
41 : * into a single vector. All overridden virtual functions are documented in
42 : * function_base.h.
43 : *
44 : * \author Roy Stogner
45 : * \date 2012
46 : * \brief Function which is a function of another function.
47 : */
48 : template <typename Output=Number>
49 : class CompositeFunction : public FunctionBase<Output>
50 : {
51 : public:
52 : explicit
53 8198 : CompositeFunction () = default;
54 :
55 : /**
56 : * This class can be default move constructed and assigned.
57 : */
58 : CompositeFunction (CompositeFunction &&) = default;
59 : CompositeFunction & operator= (CompositeFunction &&) = default;
60 :
61 : /**
62 : * This class contains unique_ptr members so it can't be default
63 : * copied or assigned.
64 : */
65 : CompositeFunction (const CompositeFunction &) = delete;
66 : CompositeFunction & operator= (const CompositeFunction &) = delete;
67 :
68 : /**
69 : * The subfunctions vector is automatically cleaned up.
70 : */
71 16628 : virtual ~CompositeFunction () = default;
72 :
73 : /**
74 : * Attach a new subfunction, along with a map from the indices of
75 : * the attached subfunction to the indices of the composed function.
76 : *
77 : * The composed function will return a vector whose value at index
78 : * \p index_map[i] is the value of the attached function at index i,
79 : * i.e.,
80 : * (*this)(x, t)(index_map[i]) will return f(x, t)(i).
81 : */
82 8198 : void attach_subfunction (const FunctionBase<Output> & f,
83 : std::vector<unsigned int> index_map)
84 : {
85 232 : const unsigned int subfunction_index =
86 464 : cast_int<unsigned int>(subfunctions.size());
87 232 : libmesh_assert_equal_to(subfunctions.size(), index_maps.size());
88 :
89 16164 : subfunctions.push_back(f.clone());
90 :
91 8198 : unsigned int max_index =
92 232 : *std::max_element(index_map.begin(), index_map.end());
93 :
94 8430 : if (max_index >= reverse_index_map.size())
95 : reverse_index_map.resize
96 10080314 : (max_index+1, std::make_pair(libMesh::invalid_uint,
97 : libMesh::invalid_uint));
98 :
99 23912 : for (auto j : index_range(index_map))
100 : {
101 444 : libmesh_assert_less(index_map[j], reverse_index_map.size());
102 444 : libmesh_assert_equal_to(reverse_index_map[index_map[j]].first,
103 : libMesh::invalid_uint);
104 444 : libmesh_assert_equal_to(reverse_index_map[index_map[j]].second,
105 : libMesh::invalid_uint);
106 16158 : reverse_index_map[index_map[j]] =
107 888 : std::make_pair(subfunction_index, j);
108 : }
109 :
110 : // Now check for time dependence
111 : // We only check the function we just added instead of researching all subfunctions
112 : // If this is the first subfunction, then that determines the time-dependence.
113 8430 : if (subfunctions.size() == 1)
114 8198 : this->_is_time_dependent = subfunctions[0]->is_time_dependent();
115 :
116 : // Otherwise, we have more than 1 function already.
117 : // If _is_time_dependent is true, then one of the previous
118 : // subfunctions is time-dependent and thus this CompositeFunction
119 : // time-dependent. If _is_time_dependent is false, then the subfunction
120 : // just added determines the time-dependence.
121 0 : else if (!this->_is_time_dependent)
122 0 : this->_is_time_dependent = (subfunctions.back())->is_time_dependent();
123 :
124 8198 : index_maps.push_back(std::move(index_map));
125 8198 : }
126 :
127 0 : virtual Output operator() (const Point & p,
128 : const Real time = 0) override
129 : {
130 0 : return this->component(0,p,time);
131 : }
132 :
133 0 : virtual void operator() (const Point & p,
134 : const Real time,
135 : DenseVector<Output> & output) override
136 : {
137 0 : libmesh_assert_greater_equal (output.size(),
138 : reverse_index_map.size());
139 :
140 : // Necessary in case we have output components not covered by
141 : // any subfunctions
142 0 : output.zero();
143 :
144 0 : DenseVector<Output> temp;
145 0 : for (auto i : index_range(subfunctions))
146 : {
147 0 : temp.resize(cast_int<unsigned int>(index_maps[i].size()));
148 0 : (*subfunctions[i])(p, time, temp);
149 0 : for (auto j : index_range(temp))
150 0 : output(index_maps[i][j]) = temp(j);
151 : }
152 0 : }
153 :
154 266853 : virtual Output component (unsigned int i,
155 : const Point & p,
156 : Real time) override
157 : {
158 312301 : if (i >= reverse_index_map.size() ||
159 266853 : reverse_index_map[i].first == libMesh::invalid_uint)
160 0 : return 0;
161 :
162 22724 : libmesh_assert_less(reverse_index_map[i].first,
163 : subfunctions.size());
164 22724 : libmesh_assert_not_equal_to(reverse_index_map[i].second,
165 : libMesh::invalid_uint);
166 266853 : return subfunctions[reverse_index_map[i].first]->
167 266853 : component(reverse_index_map[i].second,p,time);
168 : }
169 :
170 4099 : virtual std::unique_ptr<FunctionBase<Output>> clone() const override
171 : {
172 4215 : auto returnval = std::make_unique<CompositeFunction>();
173 8198 : for (auto i : index_range(subfunctions))
174 8198 : returnval->attach_subfunction(*subfunctions[i], index_maps[i]);
175 4215 : return returnval;
176 3867 : }
177 :
178 : unsigned int n_subfunctions () const
179 : {
180 : return subfunctions.size();
181 : }
182 :
183 : unsigned int n_components () const
184 : {
185 : return reverse_index_map.size();
186 : }
187 :
188 : private:
189 : // list of functions which fill in our values
190 : std::vector<std::unique_ptr<FunctionBase<Output>>> subfunctions;
191 :
192 : // for each function, list of which global indices it fills in
193 : std::vector<std::vector<unsigned int>> index_maps;
194 :
195 : // for each global index, which local index of which function is it?
196 : std::vector<std::pair<unsigned int, unsigned int>> reverse_index_map;
197 : };
198 :
199 :
200 : } // namespace libMesh
201 :
202 : #endif // LIBMESH_COMPOSITE_FUNCTION_H
|
