Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

719 lines (608 sloc) 26.377 kb
// The libMesh Finite Element Library.
// Copyright (C) 2002-2012 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
// C++ includes
#include <algorithm> // for std::fill
#include <sstream>
#include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
#include <cmath> // for sqrt
// Local Includes
#include "libmesh/dof_map.h"
#include "libmesh/elem.h"
#include "libmesh/equation_systems.h"
#include "libmesh/error_vector.h"
#include "libmesh/fe.h"
#include "libmesh/libmesh_common.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/mesh_base.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/quadrature.h"
#include "libmesh/system.h"
#include "libmesh/uniform_refinement_estimator.h"
#include "libmesh/partitioner.h"
#include "libmesh/tensor_tools.h"
#ifdef LIBMESH_ENABLE_AMR
namespace libMesh
{
//-----------------------------------------------------------------
// ErrorEstimator implementations
void UniformRefinementEstimator::estimate_error (const System& _system,
ErrorVector& error_per_cell,
const NumericVector<Number>* solution_vector,
bool estimate_parent_error)
{
START_LOG("estimate_error()", "UniformRefinementEstimator");
std::map<const System*, const NumericVector<Number>* > solution_vectors;
solution_vectors[&_system] = solution_vector;
this->_estimate_error (NULL, &_system, &error_per_cell, NULL, NULL,
&solution_vectors, estimate_parent_error);
STOP_LOG("estimate_error()", "UniformRefinementEstimator");
}
void UniformRefinementEstimator::estimate_errors (const EquationSystems& _es,
ErrorVector& error_per_cell,
const std::map<const System*, SystemNorm>& error_norms,
const std::map<const System*, const NumericVector<Number>* >* solution_vectors,
bool estimate_parent_error)
{
START_LOG("estimate_errors()", "UniformRefinementEstimator");
this->_estimate_error (&_es, NULL, &error_per_cell, NULL,
&error_norms, solution_vectors,
estimate_parent_error);
STOP_LOG("estimate_errors()", "UniformRefinementEstimator");
}
void UniformRefinementEstimator::estimate_errors (const EquationSystems& _es,
ErrorMap& errors_per_cell,
const std::map<const System*, const NumericVector<Number>* >* solution_vectors,
bool estimate_parent_error)
{
START_LOG("estimate_errors()", "UniformRefinementEstimator");
this->_estimate_error (&_es, NULL, NULL, &errors_per_cell, NULL,
solution_vectors, estimate_parent_error);
STOP_LOG("estimate_errors()", "UniformRefinementEstimator");
}
void UniformRefinementEstimator::_estimate_error (const EquationSystems* _es,
const System* _system,
ErrorVector* error_per_cell,
ErrorMap* errors_per_cell,
const std::map<const System*, SystemNorm > *_error_norms,
const std::map<const System*, const NumericVector<Number>* >* solution_vectors,
bool)
{
// Get a vector of the Systems we're going to work on,
// and set up a error_norms map if necessary
std::vector<System *> system_list;
AutoPtr<std::map<const System*, SystemNorm > > error_norms =
AutoPtr<std::map<const System*, SystemNorm > >
(new std::map<const System*, SystemNorm>);
if (_es)
{
libmesh_assert(!_system);
libmesh_assert(_es->n_systems());
_system = &(_es->get_system(0));
libmesh_assert_equal_to (&(_system->get_equation_systems()), _es);
libmesh_assert(_es->n_systems());
for (unsigned int i=0; i != _es->n_systems(); ++i)
// We have to break the rules here, because we can't refine a const System
system_list.push_back(const_cast<System *>(&(_es->get_system(i))));
// If we're computing one vector, we need to know how to scale
// each variable's contributions to it.
if (_error_norms)
{
libmesh_assert(!errors_per_cell);
}
else
// If we're computing many vectors, we just need to know which
// variables to skip
{
libmesh_assert (errors_per_cell);
_error_norms = error_norms.get();
for (unsigned int i=0; i!= _es->n_systems(); ++i)
{
const System &sys = _es->get_system(i);
unsigned int n_vars = sys.n_vars();
std::vector<Real> weights(n_vars, 0.0);
for (unsigned int v = 0; v != n_vars; ++v)
{
if (errors_per_cell->find(std::make_pair(&sys, v)) ==
errors_per_cell->end())
continue;
weights[v] = 1.0;
}
(*error_norms)[&sys] =
SystemNorm(std::vector<FEMNormType>(n_vars, error_norm.type(0)),
weights);
}
}
}
else
{
libmesh_assert(_system);
// We have to break the rules here, because we can't refine a const System
system_list.push_back(const_cast<System *>(_system));
libmesh_assert(!_error_norms);
(*error_norms)[_system] = error_norm;
_error_norms = error_norms.get();
}
// An EquationSystems reference will be convenient.
// We have to break the rules here, because we can't refine a const System
EquationSystems& es =
const_cast<EquationSystems &>(_system->get_equation_systems());
// The current mesh
MeshBase& mesh = es.get_mesh();
// The dimensionality of the mesh
const unsigned int dim = mesh.mesh_dimension();
// Resize the error_per_cell vectors to be
// the number of elements, initialize them to 0.
if (error_per_cell)
{
error_per_cell->clear();
error_per_cell->resize (mesh.max_elem_id(), 0.);
}
else
{
libmesh_assert(errors_per_cell);
for (ErrorMap::iterator i = errors_per_cell->begin();
i != errors_per_cell->end(); ++i)
{
ErrorVector *e = i->second;
e->clear();
e->resize(mesh.max_elem_id(), 0.);
}
}
// We'll want to back up all coarse grid vectors
std::vector<std::map<std::string, NumericVector<Number> *> >
coarse_vectors(system_list.size());
std::vector<NumericVector<Number> *>
coarse_solutions(system_list.size());
std::vector<NumericVector<Number> *>
coarse_local_solutions(system_list.size());
// And make copies of projected solutions
std::vector<NumericVector<Number> *>
projected_solutions(system_list.size());
// And we'll need to temporarily change solution projection settings
std::vector<bool> old_projection_settings(system_list.size());
// And it'll be best to avoid any repartitioning
AutoPtr<Partitioner> old_partitioner = mesh.partitioner();
mesh.partitioner().reset(NULL);
for (unsigned int i=0; i != system_list.size(); ++i)
{
System &system = *system_list[i];
// Check for valid error_norms
libmesh_assert (_error_norms->find(&system) !=
_error_norms->end());
// Back up the solution vector
coarse_solutions[i] = system.solution->clone().release();
coarse_local_solutions[i] =
system.current_local_solution->clone().release();
// Back up all other coarse grid vectors
for (System::vectors_iterator vec = system.vectors_begin(); vec !=
system.vectors_end(); ++vec)
{
// The (string) name of this vector
const std::string& var_name = vec->first;
coarse_vectors[i][var_name] = vec->second->clone().release();
}
// Use a non-standard solution vector if necessary
if (solution_vectors &&
solution_vectors->find(&system) != solution_vectors->end() &&
solution_vectors->find(&system)->second &&
solution_vectors->find(&system)->second != system.solution.get())
{
NumericVector<Number>* newsol =
const_cast<NumericVector<Number>*>
(solution_vectors->find(&system)->second);
newsol->swap(*system.solution);
system.update();
}
// Make sure the solution is projected when we refine the mesh
old_projection_settings[i] = system.project_solution_on_reinit();
system.project_solution_on_reinit() = true;
}
// Find the number of coarse mesh elements, to make it possible
// to find correct coarse elem ids later
const unsigned int max_coarse_elem_id = mesh.max_elem_id();
#ifndef NDEBUG
// n_coarse_elem is only used in an assertion later so
// avoid declaring it unless asserts are active.
const unsigned int n_coarse_elem = mesh.n_elem();
#endif
// Uniformly refine the mesh
MeshRefinement mesh_refinement(mesh);
libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
// FIXME: this may break if there is more than one System
// on this mesh but estimate_error was still called instead of
// estimate_errors
for (unsigned int i = 0; i != number_h_refinements; ++i)
{
mesh_refinement.uniformly_refine(1);
es.reinit();
}
for (unsigned int i = 0; i != number_p_refinements; ++i)
{
mesh_refinement.uniformly_p_refine(1);
es.reinit();
}
for (unsigned int i=0; i != system_list.size(); ++i)
{
System &system = *system_list[i];
// Copy the projected coarse grid solutions, which will be
// overwritten by solve()
// projected_solutions[i] = system.solution->clone().release();
projected_solutions[i] = NumericVector<Number>::build().release();
projected_solutions[i]->init(system.solution->size(), true, SERIAL);
system.solution->localize(*projected_solutions[i],
system.get_dof_map().get_send_list());
}
// Are we doing a forward or an adjoint solve?
bool solve_adjoint = false;
if (solution_vectors)
{
System *sys = system_list[0];
libmesh_assert (solution_vectors->find(sys) !=
solution_vectors->end());
const NumericVector<Number> *vec = solution_vectors->find(sys)->second;
for (unsigned int j=0; j != sys->qoi.size(); ++j)
{
std::ostringstream adjoint_name;
adjoint_name << "adjoint_solution" << j;
if (vec == sys->request_vector(adjoint_name.str()))
{
solve_adjoint = true;
break;
}
}
}
// Get the uniformly refined solution.
if (_es)
{
// Even if we had a decent preconditioner, valid matrix etc. before
// refinement, we don't any more.
for (unsigned int i=0; i != es.n_systems(); ++i)
es.get_system(i).disable_cache();
// No specified vectors == forward solve
if (!solution_vectors)
es.solve();
else
{
libmesh_assert_equal_to (solution_vectors->size(), es.n_systems());
libmesh_assert (solution_vectors->find(system_list[0]) !=
solution_vectors->end());
libmesh_assert(solve_adjoint ||
(solution_vectors->find(system_list[0])->second ==
system_list[0]->solution.get()) ||
!solution_vectors->find(system_list[0])->second);
#ifdef DEBUG
for (unsigned int i=0; i != system_list.size(); ++i)
{
System *sys = system_list[i];
libmesh_assert (solution_vectors->find(sys) !=
solution_vectors->end());
const NumericVector<Number> *vec = solution_vectors->find(sys)->second;
if (solve_adjoint)
{
bool found_vec = false;
for (unsigned int j=0; j != sys->qoi.size(); ++j)
{
std::ostringstream adjoint_name;
adjoint_name << "adjoint_solution" << j;
if (vec == sys->request_vector(adjoint_name.str()))
{
found_vec = true;
break;
}
}
libmesh_assert(found_vec);
}
else
libmesh_assert(vec == sys->solution.get() || !vec);
}
#endif
if (solve_adjoint)
{
std::vector<unsigned int> adjs(system_list.size(),
libMesh::invalid_uint);
// Set up proper initial guesses
for (unsigned int i=0; i != system_list.size(); ++i)
{
System *sys = system_list[i];
libmesh_assert (solution_vectors->find(sys) !=
solution_vectors->end());
const NumericVector<Number> *vec = solution_vectors->find(sys)->second;
for (unsigned int j=0; j != sys->qoi.size(); ++j)
{
std::ostringstream adjoint_name;
adjoint_name << "adjoint_solution" << j;
if (vec == sys->request_vector(adjoint_name.str()))
{
adjs[i] = j;
break;
}
}
libmesh_assert_not_equal_to (adjs[i], libMesh::invalid_uint);
system_list[i]->get_adjoint_solution(adjs[i]) =
*system_list[i]->solution;
}
es.adjoint_solve();
// Put the adjoint_solution into solution for
// comparisons
for (unsigned int i=0; i != system_list.size(); ++i)
{
system_list[i]->get_adjoint_solution(adjs[i]).swap(*system_list[i]->solution);
system_list[i]->update();
}
}
else
es.solve();
}
}
else
{
System *sys = system_list[0];
// Even if we had a decent preconditioner, valid matrix etc. before
// refinement, we don't any more.
sys->disable_cache();
// No specified vectors == forward solve
if (!solution_vectors)
sys->solve();
else
{
libmesh_assert (solution_vectors->find(sys) !=
solution_vectors->end());
const NumericVector<Number> *vec = solution_vectors->find(sys)->second;
libmesh_assert(solve_adjoint ||
(solution_vectors->find(sys)->second ==
sys->solution.get()) ||
!solution_vectors->find(sys)->second);
if (solve_adjoint)
{
unsigned int adj = libMesh::invalid_uint;
for (unsigned int j=0; j != sys->qoi.size(); ++j)
{
std::ostringstream adjoint_name;
adjoint_name << "adjoint_solution" << j;
if (vec == sys->request_vector(adjoint_name.str()))
{
adj = j;
break;
}
}
libmesh_assert_not_equal_to (adj, libMesh::invalid_uint);
// Set up proper initial guess
sys->get_adjoint_solution(adj) = *sys->solution;
sys->adjoint_solve();
// Put the adjoint_solution into solution for
// comparisons
sys->get_adjoint_solution(adj).swap(*sys->solution);
sys->update();
}
else
sys->solve();
}
}
// Get the error in the uniformly refined solution(s).
for (unsigned int i=0; i != system_list.size(); ++i)
{
System &system = *system_list[i];
unsigned int n_vars = system.n_vars();
DofMap &dof_map = system.get_dof_map();
const SystemNorm &system_i_norm =
_error_norms->find(&system)->second;
NumericVector<Number> *projected_solution = projected_solutions[i];
// Loop over all the variables in the system
for (unsigned int var=0; var<n_vars; var++)
{
// Get the error vector to fill for this system and variable
ErrorVector *err_vec = error_per_cell;
if (!err_vec)
{
libmesh_assert(errors_per_cell);
err_vec =
(*errors_per_cell)[std::make_pair(&system,var)];
}
// The type of finite element to use for this variable
const FEType& fe_type = dof_map.variable_type (var);
// Finite element object for each fine element
AutoPtr<FEBase> fe (FEBase::build (dim, fe_type));
// Build and attach an appropriate quadrature rule
AutoPtr<QBase> qrule = fe_type.default_quadrature_rule(dim);
fe->attach_quadrature_rule (qrule.get());
const std::vector<Real>& JxW = fe->get_JxW();
const std::vector<std::vector<Real> >& phi = fe->get_phi();
const std::vector<std::vector<RealGradient> >& dphi =
fe->get_dphi();
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
const std::vector<std::vector<RealTensor> >& d2phi =
fe->get_d2phi();
#endif
// The global DOF indices for the fine element
std::vector<unsigned int> dof_indices;
// Iterate over all the active elements in the fine mesh
// that live on this processor.
MeshBase::const_element_iterator elem_it = mesh.active_local_elements_begin();
const MeshBase::const_element_iterator elem_end = mesh.active_local_elements_end();
for (; elem_it != elem_end; ++elem_it)
{
// e is necessarily an active element on the local processor
const Elem* elem = *elem_it;
// Find the element id for the corresponding coarse grid element
const Elem* coarse = elem;
unsigned int e_id = coarse->id();
while (e_id >= max_coarse_elem_id)
{
libmesh_assert (coarse->parent());
coarse = coarse->parent();
e_id = coarse->id();
}
double L2normsq = 0., H1seminormsq = 0., H2seminormsq = 0.;
// reinitialize the element-specific data
// for the current element
fe->reinit (elem);
// Get the local to global degree of freedom maps
dof_map.dof_indices (elem, dof_indices, var);
// The number of quadrature points
const unsigned int n_qp = qrule->n_points();
// The number of shape functions
const unsigned int n_sf = dof_indices.size();
//
// Begin the loop over the Quadrature points.
//
for (unsigned int qp=0; qp<n_qp; qp++)
{
Number u_fine = 0., u_coarse = 0.;
Gradient grad_u_fine, grad_u_coarse;
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
Tensor grad2_u_fine, grad2_u_coarse;
#endif
// Compute solution values at the current
// quadrature point. This reqiures a sum
// over all the shape functions evaluated
// at the quadrature point.
for (unsigned int i=0; i<n_sf; i++)
{
u_fine += phi[i][qp]*system.current_solution (dof_indices[i]);
u_coarse += phi[i][qp]*(*projected_solution) (dof_indices[i]);
grad_u_fine += dphi[i][qp]*system.current_solution (dof_indices[i]);
grad_u_coarse += dphi[i][qp]*(*projected_solution) (dof_indices[i]);
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
grad2_u_fine += d2phi[i][qp]*system.current_solution (dof_indices[i]);
grad2_u_coarse += d2phi[i][qp]*(*projected_solution) (dof_indices[i]);
#endif
}
// Compute the value of the error at this quadrature point
const Number val_error = u_fine - u_coarse;
// Add the squares of the error to each contribution
if (system_i_norm.type(var) == L2 ||
system_i_norm.type(var) == H1 ||
system_i_norm.type(var) == H2)
{
L2normsq += JxW[qp] * system_i_norm.weight_sq(var) *
TensorTools::norm_sq(val_error);
libmesh_assert_greater_equal (L2normsq, 0.);
}
// Compute the value of the error in the gradient at this
// quadrature point
if (system_i_norm.type(var) == H1 ||
system_i_norm.type(var) == H2 ||
system_i_norm.type(var) == H1_SEMINORM)
{
Gradient grad_error = grad_u_fine - grad_u_coarse;
H1seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
grad_error.size_sq();
libmesh_assert_greater_equal (H1seminormsq, 0.);
}
#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
// Compute the value of the error in the hessian at this
// quadrature point
if (system_i_norm.type(var) == H2 ||
system_i_norm.type(var) == H2_SEMINORM)
{
Tensor grad2_error = grad2_u_fine - grad2_u_coarse;
H2seminormsq += JxW[qp] * system_i_norm.weight_sq(var) *
grad2_error.size_sq();
libmesh_assert_greater_equal (H2seminormsq, 0.);
}
#endif
} // end qp loop
if (system_i_norm.type(var) == L2 ||
system_i_norm.type(var) == H1 ||
system_i_norm.type(var) == H2)
(*err_vec)[e_id] += L2normsq;
if (system_i_norm.type(var) == H1 ||
system_i_norm.type(var) == H2 ||
system_i_norm.type(var) == H1_SEMINORM)
(*err_vec)[e_id] += H1seminormsq;
if (system_i_norm.type(var) == H2 ||
system_i_norm.type(var) == H2_SEMINORM)
(*err_vec)[e_id] += H2seminormsq;
} // End loop over active local elements
} // End loop over variables
// Don't bother projecting the solution; we'll restore from backup
// after coarsening
system.project_solution_on_reinit() = false;
}
// Uniformly coarsen the mesh, without projecting the solution
libmesh_assert (number_h_refinements > 0 || number_p_refinements > 0);
for (unsigned int i = 0; i != number_h_refinements; ++i)
{
mesh_refinement.uniformly_coarsen(1);
// FIXME - should the reinits here be necessary? - RHS
es.reinit();
}
for (unsigned int i = 0; i != number_p_refinements; ++i)
{
mesh_refinement.uniformly_p_coarsen(1);
es.reinit();
}
// We should be back where we started
libmesh_assert_equal_to (n_coarse_elem, mesh.n_elem());
// Each processor has now computed the error contribuions
// for its local elements. We need to sum the vector
// and then take the square-root of each component. Note
// that we only need to sum if we are running on multiple
// processors, and we only need to take the square-root
// if the value is nonzero. There will in general be many
// zeros for the inactive elements.
if (error_per_cell)
{
// First sum the vector of estimated error values
this->reduce_error(*error_per_cell);
// Compute the square-root of each component.
START_LOG("std::sqrt()", "UniformRefinementEstimator");
for (unsigned int i=0; i<error_per_cell->size(); i++)
if ((*error_per_cell)[i] != 0.)
(*error_per_cell)[i] = std::sqrt((*error_per_cell)[i]);
STOP_LOG("std::sqrt()", "UniformRefinementEstimator");
}
else
{
for (ErrorMap::iterator i = errors_per_cell->begin();
i != errors_per_cell->end(); ++i)
{
ErrorVector *e = i->second;
// First sum the vector of estimated error values
this->reduce_error(*e);
// Compute the square-root of each component.
START_LOG("std::sqrt()", "UniformRefinementEstimator");
for (unsigned int i=0; i<e->size(); i++)
if ((*e)[i] != 0.)
(*e)[i] = std::sqrt((*e)[i]);
STOP_LOG("std::sqrt()", "UniformRefinementEstimator");
}
}
// Restore old solutions and clean up the heap
for (unsigned int i=0; i != system_list.size(); ++i)
{
System &system = *system_list[i];
system.project_solution_on_reinit() = old_projection_settings[i];
// Restore the coarse solution vectors and delete their copies
*system.solution = *coarse_solutions[i];
delete coarse_solutions[i];
*system.current_local_solution = *coarse_local_solutions[i];
delete coarse_local_solutions[i];
delete projected_solutions[i];
for (System::vectors_iterator vec = system.vectors_begin(); vec !=
system.vectors_end(); ++vec)
{
// The (string) name of this vector
const std::string& var_name = vec->first;
system.get_vector(var_name) = *coarse_vectors[i][var_name];
coarse_vectors[i][var_name]->clear();
delete coarse_vectors[i][var_name];
}
}
// Restore old partitioner settings
mesh.partitioner() = old_partitioner;
}
} // namespace libMesh
#endif // #ifdef LIBMESH_ENABLE_AMR
Jump to Line
Something went wrong with that request. Please try again.