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             : 
      19             : // Local include files
      20             : #include "libmesh/libmesh_common.h"
      21             : #include "libmesh/error_estimator.h"
      22             : #include "libmesh/error_vector.h"
      23             : #include "libmesh/equation_systems.h"
      24             : #include "libmesh/parallel.h"
      25             : #include "libmesh/enum_error_estimator_type.h"
      26             : #include "libmesh/int_range.h"
      27             : 
      28             : namespace libMesh
      29             : {
      30             : 
      31             : // ErrorEstimator functions
      32       52728 : void ErrorEstimator::reduce_error (std::vector<ErrorVectorReal> & error_per_cell,
      33             :                                    const Parallel::Communicator & comm) const
      34             : {
      35             :   // This function must be run on all processors at once
      36             :   // parallel_object_only();
      37             : 
      38             :   // Each processor has now computed the error contributions
      39             :   // for its local elements.  We may need to sum the vector to
      40             :   // recover the error for each element.
      41             : 
      42       52728 :   comm.sum(error_per_cell);
      43       52728 : }
      44             : 
      45             : 
      46             : 
      47           0 : void ErrorEstimator::estimate_errors(const EquationSystems & equation_systems,
      48             :                                      ErrorVector & error_per_cell,
      49             :                                      const std::map<const System *, SystemNorm> & error_norms,
      50             :                                      const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
      51             :                                      bool estimate_parent_error)
      52             : {
      53           0 :   SystemNorm old_error_norm = this->error_norm;
      54             : 
      55             :   // Sum the error values from each system
      56           0 :   for (auto s : make_range(equation_systems.n_systems()))
      57             :     {
      58           0 :       ErrorVector system_error_per_cell;
      59           0 :       const System & sys = equation_systems.get_system(s);
      60           0 :       if (const auto it = error_norms.find(&sys);
      61           0 :           it == error_norms.end())
      62           0 :         this->error_norm = old_error_norm;
      63             :       else
      64           0 :         this->error_norm = it->second;
      65             : 
      66           0 :       const NumericVector<Number> * solution_vector = nullptr;
      67           0 :       if (solution_vectors &&
      68           0 :           solution_vectors->find(&sys) != solution_vectors->end())
      69           0 :         solution_vector = solution_vectors->find(&sys)->second;
      70             : 
      71           0 :       this->estimate_error(sys, system_error_per_cell,
      72           0 :                            solution_vector, estimate_parent_error);
      73             : 
      74           0 :       if (s)
      75             :         {
      76           0 :           libmesh_assert_equal_to (error_per_cell.size(), system_error_per_cell.size());
      77           0 :           for (auto i : index_range(error_per_cell))
      78           0 :             error_per_cell[i] += system_error_per_cell[i];
      79             :         }
      80             :       else
      81           0 :         error_per_cell = system_error_per_cell;
      82             :     }
      83             : 
      84             :   // Restore our old state before returning
      85           0 :   this->error_norm = old_error_norm;
      86           0 : }
      87             : 
      88             : 
      89             : 
      90             : /**
      91             :  * FIXME: This is a default implementation - derived classes should
      92             :  * reimplement it for efficiency.
      93             :  */
      94        2800 : void ErrorEstimator::estimate_errors(const EquationSystems & equation_systems,
      95             :                                      ErrorMap & errors_per_cell,
      96             :                                      const std::map<const System *, const NumericVector<Number> *> * solution_vectors,
      97             :                                      bool estimate_parent_error)
      98             : {
      99        2960 :   SystemNorm old_error_norm = this->error_norm;
     100             : 
     101             :   // Find the requested error values from each system
     102        5600 :   for (auto s : make_range(equation_systems.n_systems()))
     103             :     {
     104          80 :       const System & sys = equation_systems.get_system(s);
     105             : 
     106          80 :       unsigned int n_vars = sys.n_vars();
     107             : 
     108       14000 :       for (unsigned int v = 0; v != n_vars; ++v)
     109             :         {
     110             :           // Only fill in ErrorVectors the user asks for
     111       11200 :           if (!errors_per_cell.count(std::make_pair(&sys, v)))
     112           0 :             continue;
     113             : 
     114             :           // Calculate error in only one variable
     115       11200 :           std::vector<Real> weights(n_vars, 0.0);
     116       11200 :           weights[v] = 1.0;
     117             :           this->error_norm =
     118       22400 :             SystemNorm(std::vector<FEMNormType>(n_vars, old_error_norm.type(v)),
     119       11200 :                        weights);
     120             : 
     121         320 :           const NumericVector<Number> * solution_vector = nullptr;
     122       11840 :           if (solution_vectors &&
     123       11520 :               solution_vectors->find(&sys) != solution_vectors->end())
     124       11200 :             solution_vector = solution_vectors->find(&sys)->second;
     125             : 
     126             :           this->estimate_error
     127       11200 :             (sys, *errors_per_cell[std::make_pair(&sys, v)],
     128         640 :              solution_vector, estimate_parent_error);
     129             :         }
     130             :     }
     131             : 
     132             :   // Restore our old state before returning
     133        2800 :   this->error_norm = old_error_norm;
     134        2800 : }
     135             : 
     136             : } // namespace libMesh