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_FEM_FUNCTION_H
      19             : #define LIBMESH_COMPOSITE_FEM_FUNCTION_H
      20             : 
      21             : // libMesh includes
      22             : #include "libmesh/dense_vector.h"
      23             : #include "libmesh/fem_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             :  * FEMFunction which is a function of another function.
      38             :  *
      39             :  * \author Roy Stogner
      40             :  * \date 2012
      41             :  * \brief FEMFunction which is a function of another function.
      42             :  */
      43             : template <typename Output=Number>
      44             : class CompositeFEMFunction : public FEMFunctionBase<Output>
      45             : {
      46             : public:
      47             :   explicit
      48           0 :   CompositeFEMFunction () = default;
      49             : 
      50             :   /**
      51             :    * This class can be default move constructed and assigned.
      52             :    */
      53             :   CompositeFEMFunction (CompositeFEMFunction &&) = default;
      54             :   CompositeFEMFunction & operator= (CompositeFEMFunction &&) = default;
      55             : 
      56             :   /**
      57             :    * This class contains unique_ptr members so it can't be default
      58             :    * copied or assigned.
      59             :    */
      60             :   CompositeFEMFunction (const CompositeFEMFunction &) = delete;
      61             :   CompositeFEMFunction & operator= (const CompositeFEMFunction &) = delete;
      62             : 
      63             :   /**
      64             :    * The subfunctions vector is automatically cleaned up.
      65             :    */
      66           0 :   virtual ~CompositeFEMFunction () = default;
      67             : 
      68             :   /**
      69             :    * Attach a new subfunction, along with a map from the indices of
      70             :    * that subfunction to the indices of the global function.
      71             :    * (*this)(index_map[i]) will return f(i).
      72             :    */
      73           0 :   void attach_subfunction (const FEMFunctionBase<Output> & f,
      74             :                            std::vector<unsigned int> index_map)
      75             :   {
      76           0 :     const unsigned int subfunction_index =
      77           0 :       cast_int<unsigned int>(subfunctions.size());
      78           0 :     libmesh_assert_equal_to(subfunctions.size(), index_maps.size());
      79             : 
      80           0 :     subfunctions.push_back(f.clone());
      81             : 
      82           0 :     unsigned int max_index =
      83           0 :       *std::max_element(index_map.begin(), index_map.end());
      84             : 
      85           0 :     if (max_index >= reverse_index_map.size())
      86             :       reverse_index_map.resize
      87           0 :         (max_index+1, std::make_pair(libMesh::invalid_uint,
      88             :                                      libMesh::invalid_uint));
      89             : 
      90           0 :     for (auto j : index_range(index_map))
      91             :       {
      92           0 :         libmesh_assert_less(index_map[j], reverse_index_map.size());
      93           0 :         libmesh_assert_equal_to(reverse_index_map[index_map[j]].first,
      94             :                                 libMesh::invalid_uint);
      95           0 :         libmesh_assert_equal_to(reverse_index_map[index_map[j]].second,
      96             :                                 libMesh::invalid_uint);
      97           0 :         reverse_index_map[index_map[j]] =
      98           0 :           std::make_pair(subfunction_index, j);
      99             :       }
     100             : 
     101           0 :     index_maps.push_back(std::move(index_map));
     102           0 :   }
     103             : 
     104           0 :   virtual Output operator() (const FEMContext & c,
     105             :                              const Point & p,
     106             :                              const Real time = 0) override
     107             :   {
     108           0 :     return this->component(c,0,p,time);
     109             :   }
     110             : 
     111           0 :   virtual void operator() (const FEMContext & c,
     112             :                            const Point & p,
     113             :                            const Real time,
     114             :                            DenseVector<Output> & output) override
     115             :   {
     116           0 :     libmesh_assert_greater_equal (output.size(),
     117             :                                   reverse_index_map.size());
     118             : 
     119             :     // Necessary in case we have output components not covered by
     120             :     // any subfunctions
     121           0 :     output.zero();
     122             : 
     123           0 :     DenseVector<Output> temp;
     124           0 :     for (auto i : index_range(subfunctions))
     125             :       {
     126           0 :         temp.resize(cast_int<unsigned int>(index_maps[i].size()));
     127           0 :         (*subfunctions[i])(c, p, time, temp);
     128           0 :         for (auto j : index_range(temp))
     129           0 :           output(index_maps[i][j]) = temp(j);
     130             :       }
     131           0 :   }
     132             : 
     133           0 :   virtual Output component (const FEMContext & c,
     134             :                             unsigned int i,
     135             :                             const Point & p,
     136             :                             Real time) override
     137             :   {
     138           0 :     if (i >= reverse_index_map.size() ||
     139           0 :         reverse_index_map[i].first == libMesh::invalid_uint)
     140           0 :       return 0;
     141             : 
     142           0 :     libmesh_assert_less(reverse_index_map[i].first,
     143             :                         subfunctions.size());
     144           0 :     libmesh_assert_not_equal_to(reverse_index_map[i].second,
     145             :                                 libMesh::invalid_uint);
     146           0 :     return subfunctions[reverse_index_map[i].first]->
     147           0 :       component(c, reverse_index_map[i].second, p, time);
     148             :   }
     149             : 
     150           0 :   virtual std::unique_ptr<FEMFunctionBase<Output>> clone() const override
     151             :   {
     152           0 :     auto returnval = std::make_unique<CompositeFEMFunction>();
     153           0 :     for (auto i : index_range(subfunctions))
     154           0 :       returnval->attach_subfunction(*subfunctions[i], index_maps[i]);
     155           0 :     return returnval;
     156           0 :   }
     157             : 
     158             :   unsigned int n_subfunctions () const
     159             :   {
     160             :     return subfunctions.size();
     161             :   }
     162             : 
     163             :   unsigned int n_components () const
     164             :   {
     165             :     return reverse_index_map.size();
     166             :   }
     167             : 
     168             : private:
     169             :   // list of functions which fill in our values
     170             :   std::vector<std::unique_ptr<FEMFunctionBase<Output>>> subfunctions;
     171             : 
     172             :   // for each function, list of which global indices it fills in
     173             :   std::vector<std::vector<unsigned int>> index_maps;
     174             : 
     175             :   // for each global index, which local index of which function is it?
     176             :   std::vector<std::pair<unsigned int, unsigned int>> reverse_index_map;
     177             : };
     178             : 
     179             : 
     180             : } // namespace libMesh
     181             : 
     182             : #endif // LIBMESH_COMPOSITE_FEM_FUNCTION_H