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             : // Local includes
      19             : #include "libmesh/inter_mesh_projection.h"
      20             : 
      21             : namespace libMesh
      22             : {
      23           0 : GradientMeshFunction::GradientMeshFunction(const MeshFunction & _mesh_function):
      24           0 :   mesh_function(std::make_unique<MeshFunction>(_mesh_function))
      25             : {
      26             :   libmesh_experimental();
      27           0 : }
      28             : 
      29           0 : void GradientMeshFunction::operator() (const Point & p, const Real, DenseVector<Gradient> & output)
      30             : {
      31           0 :   Real time = 0.0;
      32           0 :   mesh_function->gradient(p, time, output.get_values());
      33           0 :   return;
      34             : }
      35             : 
      36           0 : InterMeshProjection::InterMeshProjection(System & _from_system, System & _to_system) :
      37           0 :   from_system(_from_system),
      38           0 :   to_system(_to_system)
      39             : {
      40             :   libmesh_experimental();
      41           0 : }
      42             : 
      43           0 : void InterMeshProjection::project_system_vectors()
      44             : {
      45             :   // Number of vectors to be projected
      46           0 :   libmesh_assert_equal_to (to_system.n_vectors(), from_system.n_vectors());
      47             : 
      48             :   // Number of variable components in each vector
      49           0 :   unsigned int n_vars = from_system.n_vars();
      50           0 :   libmesh_assert_equal_to (to_system.n_vars(), n_vars);
      51             : 
      52             :   // We are going to use the multi-variable MeshFunction, so we can pass
      53             :   // a single vector of variables rather than have a MeshFunction for each variable
      54           0 :   std::vector<unsigned int> variables_vector;
      55             : 
      56           0 :   for (unsigned int j = 0; j != n_vars; ++j)
      57             :     {
      58           0 :       libmesh_assert_equal_to (to_system.variable_name(j), from_system.variable_name(j));
      59           0 :       variables_vector.push_back(j);
      60             :     }
      61             : 
      62             :   // Any system holds the solution along with the other vectors system.vectors
      63             :   // We will first project the solution and then move to the system.vectors
      64             : 
      65             :   // Construct local version of the current system vector
      66             :   // This has to be a serial vector
      67             :   // Roy's FIXME: Technically it just has to be a ghosted vector whose algebraically
      68             :   // ghosted values cover a domain which is a superset of the to-system's domain ...
      69             :   // that's hard to do and we can skip it until the poor scalability bites someone.
      70           0 :   std::unique_ptr<NumericVector<Number>> solution_vector_serial = NumericVector<Number>::build(from_system.comm());
      71           0 :   solution_vector_serial->init(from_system.solution->size(), true, SERIAL);
      72             : 
      73           0 :   std::vector<Number> solution_vector;
      74           0 :   from_system.update_global_solution(solution_vector);
      75           0 :   (*solution_vector_serial) = solution_vector;
      76             : 
      77             :   // Construct a MeshFunction for the solution
      78           0 :   MeshFunction mesh_func_solution(from_system.get_equation_systems(), *solution_vector_serial, from_system.get_dof_map(), variables_vector);
      79             : 
      80           0 :   mesh_func_solution.init();
      81             : 
      82             :   // For some element types (say C1) we also need to pass a gradient evaluation MeshFunction
      83             :   // To do this evaluate, a new shim class GradientMeshFunction has been added which redirects
      84             :   // gptr::operator evaluations inside projection methods into MeshFunction::gradient calls.
      85           0 :   GradientMeshFunction gptr_solution(mesh_func_solution);
      86           0 :   gptr_solution.init();
      87             : 
      88           0 :   to_system.project_vector(*to_system.solution, &mesh_func_solution, &gptr_solution);
      89             : 
      90             :   // Now loop over the vectors in system.vectors (includes old_nonlin_sol, rhs, adjoints, adjoint_rhs, sensitivity_rhs)
      91           0 :   for (System::vectors_iterator vec = from_system.vectors_begin(), vec_end = from_system.vectors_end(); vec != vec_end; ++vec)
      92             :     {
      93             :       // The name of this vector
      94           0 :       const std::string & vec_name = vec->first;
      95             : 
      96             :       // Construct local version of the current system vector
      97             :       // This has to be a serial vector
      98             :       // Roy's FIXME: Technically it just has to be a ghosted vector whose algebraically
      99             :       // ghosted values cover a domain which is a superset of the to-system's domain ...
     100             :       // that's hard to do and we can skip it until the poor scalability bites someone.
     101           0 :       std::unique_ptr<NumericVector<Number>> current_vector_proxy = NumericVector<Number>::build(from_system.comm());
     102           0 :       current_vector_proxy->init(from_system.get_vector(vec_name).size(), true, SERIAL);
     103             : 
     104           0 :       from_system.get_vector(vec_name).localize(*current_vector_proxy);
     105             : 
     106             :       // Construct a MeshFunction for the current component
     107           0 :       MeshFunction mesh_func(from_system.get_equation_systems(), *current_vector_proxy, from_system.get_dof_map(), variables_vector);
     108           0 :       mesh_func.init();
     109             : 
     110             :       // Project the current system vector, you need to check if this vector is an adjoint to pass
     111             :       // the right options to project_vector
     112           0 :       GradientMeshFunction gptr(mesh_func);
     113           0 :       gptr.init();
     114             : 
     115             :       // The fourth argument here is whether the vector is an adjoint solution or not, we will be getting that information
     116             :       // via the from_system instead of the to_system in case the user has not set the System::_vector_is_adjoint map to true.
     117           0 :       to_system.project_vector(to_system.get_vector(vec_name), &mesh_func, &gptr, from_system.vector_is_adjoint(vec_name));
     118             : 
     119           0 :     }
     120             :   // End loop over the vectors in the system
     121             : 
     122           0 : }
     123             : // End InterMeshProjection::project_system_vectors
     124             : }
     125             : // End namespace libMesh