From: Wolfgang Bangerth Date: Tue, 19 Jul 2022 21:57:45 +0000 (-0600) Subject: Add a new code gallery program. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=edb3dfd510031ab9ee4f013cae871617ea8a48ed;p=code-gallery.git Add a new code gallery program. --- diff --git a/information_based_mesh_refinement/CMakeLists.txt b/information_based_mesh_refinement/CMakeLists.txt new file mode 100644 index 0000000..ca80f3f --- /dev/null +++ b/information_based_mesh_refinement/CMakeLists.txt @@ -0,0 +1,39 @@ +## +# CMake script for the helmholtz tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "mesh_refinement") + +# Declare all source files the target consists of. Here, this is only +# the one step-X.cc file, but as you expand your project you may wish +# to add other source files as well. If your project becomes much larger, +# you may want to either replace the following statement by something like +# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc") +# FILE(GLOB_RECURSE TARGET_INC "include/*.h") +# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) +# or switch altogether to the large project CMakeLists.txt file discussed +# in the "CMake in user projects" page accessible from the "User info" +# page of the documentation. +SET(TARGET_SRC + ${TARGET}.cc + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 9.3.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/information_based_mesh_refinement/Readme.md b/information_based_mesh_refinement/Readme.md new file mode 100644 index 0000000..a4c0449 --- /dev/null +++ b/information_based_mesh_refinement/Readme.md @@ -0,0 +1,71 @@ +Readme file for "Information density-based mesh refinement" +=========================================================== + +@note This program implements the ideas and algorithms described in + the paper "Estimating and using information in inverse problems" by + Wolfgang Bangerth, Chris R. Johnson, Dennis K. Njeru, and Bart van + Bloemen Waanders, 2022. See there for more information. + +Motivation +---------- + +Inverse problems are problems where we would like to infer properties +of a system from measurements of the system's state or response to +external stimuli. The specific example this program addresses is that +we want to identify the source term (i.e., right hand side function) +in an advection-diffusion equation from point measurements of the +solution of the equation. A typical application is that we would like +to find out the locations and strengths of pollution sources based on +measuring the concentration of the polluting substance at a number of +points. + +It is clear that in order to solve such problems, one needs to "know" +something about the system's state (here: the pollution concentration) +through measurements. Intuitively, it is also clear that we know +"more" about the pollution sources by (i) measuring at more points, +and (ii) by measuring *downstream* from the sources than we would if +we had measured *upstream*. Intuitive concepts such as this motivate +wondering whether we can define an "information density" function +whose value at a point $\mathbf x$ describes how much we know about potential +sources located at $\mathbf x$. + +The paper which this code accompanies explores the concept of +information in inverse problems. It defines an "information density" +by solving auxiliary problems for each measurement, and then outlines +possible applications for these information densities in three +vignettes: spatially variable regularization; mesh refinement; and +optimal experimental design. It then considers one of these in detail +through numerical experiments, namely mesh refinement. This program +implements the algorithms shown there and produces the numerical +results. + + +To run the code +--------------- + +After running `cmake` and compiling via `make` (or, if you have used +the `-G ...` option of `cmake`, compiling the program via your +favorite integrated development environment), you can run the +executable by either just saying `make run` or using `./mesh_refinement` +on the command line. The default is to compile in "debug mode"; you +can switch to "release mode" by saying `make release` and then +compiling everything again. + +The program contains a switch that decides which mesh refinement +algorithm to use. By default, it refines the mesh based on the +information criterion discussed in the paper; it runs a sequence +of 7 mesh refinement cycles. In debug mode, running the program as +is takes about 50 CPU minutes on a reasonably modern laptop. (The +program takes about five and a half minutes in release mode.) It +parallelizes certain operations, so the actual run time may be shorter +depending on how many cores are available. + +For each cycle, it outputs the solution as a VTU file, along with the +$A$, $B$, $C$, and $M$ matrices discussed in the paper. These matrices +can then be used to compute the eigenvalues of the $H$ matrix defined +by $H = B^T A^{-T} C A^{-1} B + \beta M$ where $\beta$ is the +regularization parameters. + +Some of the pictures shown in the paper are also reproduced as part of +this code gallery program. See the paper for captions and more information. + diff --git a/information_based_mesh_refinement/doc/adjoint.png b/information_based_mesh_refinement/doc/adjoint.png new file mode 100644 index 0000000..41e8220 Binary files /dev/null and b/information_based_mesh_refinement/doc/adjoint.png differ diff --git a/information_based_mesh_refinement/doc/author b/information_based_mesh_refinement/doc/author new file mode 100644 index 0000000..a73544d --- /dev/null +++ b/information_based_mesh_refinement/doc/author @@ -0,0 +1 @@ +Wolfgang Bangerth diff --git a/information_based_mesh_refinement/doc/builds-on b/information_based_mesh_refinement/doc/builds-on new file mode 100644 index 0000000..1740273 --- /dev/null +++ b/information_based_mesh_refinement/doc/builds-on @@ -0,0 +1 @@ +step-6 diff --git a/information_based_mesh_refinement/doc/dependencies b/information_based_mesh_refinement/doc/dependencies new file mode 100644 index 0000000..c873efc --- /dev/null +++ b/information_based_mesh_refinement/doc/dependencies @@ -0,0 +1 @@ +UMFPACK diff --git a/information_based_mesh_refinement/doc/entry-name b/information_based_mesh_refinement/doc/entry-name new file mode 100644 index 0000000..eb32993 --- /dev/null +++ b/information_based_mesh_refinement/doc/entry-name @@ -0,0 +1 @@ +Information density-based mesh refinement diff --git a/information_based_mesh_refinement/doc/recovered-forward-solution.png b/information_based_mesh_refinement/doc/recovered-forward-solution.png new file mode 100644 index 0000000..1896403 Binary files /dev/null and b/information_based_mesh_refinement/doc/recovered-forward-solution.png differ diff --git a/information_based_mesh_refinement/doc/recovered-parameter.png b/information_based_mesh_refinement/doc/recovered-parameter.png new file mode 100644 index 0000000..099129c Binary files /dev/null and b/information_based_mesh_refinement/doc/recovered-parameter.png differ diff --git a/information_based_mesh_refinement/doc/synthetic.png b/information_based_mesh_refinement/doc/synthetic.png new file mode 100644 index 0000000..96fd73f Binary files /dev/null and b/information_based_mesh_refinement/doc/synthetic.png differ diff --git a/information_based_mesh_refinement/doc/tooltip b/information_based_mesh_refinement/doc/tooltip new file mode 100644 index 0000000..f421b0a --- /dev/null +++ b/information_based_mesh_refinement/doc/tooltip @@ -0,0 +1 @@ +Implementation of a mesh refinement criterion for an inverse source identification problem based on an "information density". diff --git a/information_based_mesh_refinement/mesh_refinement.cc b/information_based_mesh_refinement/mesh_refinement.cc new file mode 100644 index 0000000..895c4c9 --- /dev/null +++ b/information_based_mesh_refinement/mesh_refinement.cc @@ -0,0 +1,889 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2022 by the deal.II authors and Wolfgang Bangerth. + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, 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. + * The full text of the license can be found in the file LICENSE.md at + * the top level directory of deal.II. + * + * --------------------------------------------------------------------- + + * + * Author: Wolfgang Bangerth, Colorado State University, 2022. + */ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +using namespace dealii; + + +// The following is the main class. It resembles a variation of the step-6 +// principal class, with the addition of information-specific stuff. It also +// has to deal with solving a vector-valued problem for (c,lambda,f) as +// primal variable, dual variable, and right hand side, as explained +// in the paper. +template +class InformationDensityMeshRefinement +{ +public: + InformationDensityMeshRefinement (); + void run (); + +private: + void compute_synthetic_measurements(); + void bounce_measurement_points_to_cell_centers (); + void setup_system(); + void assemble_system (); + void solve (); + void compute_information_content (); + void output_results (const unsigned int cycle) const; + void refine_grid (); + + const Point source_location; + const double source_radius; + + std::vector> detector_locations; + + const double regularization_parameter; + Tensor<1,dim> velocity; + + Triangulation triangulation; + FESystem fe; + DoFHandler dof_handler; + + AffineConstraints hanging_node_constraints; + + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix system_matrix; + + BlockVector solution; + BlockVector system_rhs; + + Vector information_content; + + std::vector> detector_locations_on_mesh; + std::vector measurement_values; + std::vector noise_level; +}; + + + + +template +InformationDensityMeshRefinement::InformationDensityMeshRefinement () +: +source_location (Point(-0.25,0)), +source_radius (0.2), +regularization_parameter (10000), +fe (FE_Q(3), 1, // c + FE_Q(3), 1, // lambda + FE_DGQ(0), 1), // f +dof_handler (triangulation) +{ + velocity[0] = 100; + + // We have 50 detector points on an outer ring... + for (unsigned int i=0; i<50; ++i) + { + const Point p (0.6 * std::sin(2*numbers::PI * i/50), + 0.6 * std::cos(2*numbers::PI * i/50)); + detector_locations.push_back (p); + } + + // ...and another 50 detector points on an innner ring: + for (unsigned int i=0; i<50; ++i) + { + const Point p (0.2 * std::sin(2*numbers::PI * i/50), + 0.2 * std::cos(2*numbers::PI * i/50)); + detector_locations.push_back (p); + } + + // Generate the grid we will work on: + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (4); + + // The detector locations are static, so we can already here + // generate a file that contains their locations. We use the + // particle framework to do this, using detector locations as + // particle locations. + { + Particles::ParticleHandler particle_handler(triangulation, + StaticMappingQ1::mapping); + for (const auto &loc : detector_locations) + { + Particles::Particle new_particle; + new_particle.set_location(loc); + // Insert the particle. It is a lie that the particle is in + // the first cell, but nothing we do actually cares about the + // cell a particle is in. + particle_handler.insert_particle(new_particle, + triangulation.begin_active()); + } + + Particles::DataOut particle_out; + particle_out.build_patches(particle_handler); + std::ofstream output("detector_locations.vtu"); + particle_out.write_vtu(output); + } + + // While we're generating output, also output the source location. Do this + // by outputting many (1000) points that indicate the perimeter of the source + { + Particles::ParticleHandler particle_handler(triangulation, + StaticMappingQ1::mapping); + + const unsigned int n_points = 1000; + for (unsigned int i=0; i loc = source_location; + loc[0] += source_radius * std::cos(2*numbers::PI*i/n_points); + loc[1] += source_radius * std::sin(2*numbers::PI*i/n_points); + + Particles::Particle new_particle; + new_particle.set_location(loc); + particle_handler.insert_particle(new_particle, + triangulation.begin_active()); + } + + Particles::DataOut particle_out; + particle_out.build_patches(particle_handler); + std::ofstream output("source_locations.vtu"); + particle_out.write_vtu(output); + } +} + + + +// The following function solves a forward problem on a twice +// refined mesh to compute "synthetic data". Refining the mesh +// beyond the mesh used for the inverse problem avoids an +// inverse crime. +template +void InformationDensityMeshRefinement::compute_synthetic_measurements () +{ + std::cout << "Computing synthetic data by solving the forward problem..." + << std::flush; + + // Create a triangulation and DoFHandler that corresponds to a + // twice-refined mesh so that we obtain the synthetic data with + // higher accuracy than we do on the regular mesh used for all other + // computations. + Triangulation forward_triangulation; + forward_triangulation.copy_triangulation (triangulation); + forward_triangulation.refine_global (2); + + const FE_Q forward_fe (fe.base_element(0).degree); + DoFHandler forward_dof_handler (forward_triangulation); + forward_dof_handler.distribute_dofs (forward_fe); + + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(forward_dof_handler, constraints); + constraints.close(); + + SparsityPattern sparsity (forward_dof_handler.n_dofs(), + forward_dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (forward_dof_handler, sparsity); + constraints.condense (sparsity); + sparsity.compress (); + + SparseMatrix system_matrix (sparsity); + Vector system_rhs (forward_dof_handler.n_dofs()); + + QGauss quadrature_formula(3); + FEValues fe_values (forward_fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = forward_fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + // First assemble the system matrix and right hand side for the forward + // problem: + { + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + for (const auto &cell : forward_dof_handler.active_cell_iterators()) + { + fe_values.reinit (cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + constraints.distribute_local_to_global (cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + + std::map boundary_values; + VectorTools::interpolate_boundary_values (forward_dof_handler, + 0, + Functions::ZeroFunction(), + boundary_values); + Vector tmp (forward_dof_handler.n_dofs()); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + tmp, + system_rhs); + } + + // Solve the forward problem and output it into its own VTU file: + SparseDirectUMFPACK A_inverse; + Vector forward_solution (forward_dof_handler.n_dofs()); + forward_solution = system_rhs; + A_inverse.solve(system_matrix, forward_solution); + + const double max_forward_solution = forward_solution.linfty_norm(); + + { + DataOut data_out; + data_out.attach_dof_handler (forward_dof_handler); + data_out.add_data_vector (forward_solution, "c"); + data_out.build_patches (4); + + std::ofstream out ("forward-solution.vtu"); + data_out.write_vtu (out); + } + + // Now evaluate the forward solution at the measurement points: + for (const auto &p : detector_locations) + { + // same 10% noise level for all points + noise_level.push_back (0.1 * max_forward_solution); + + const double z_n = VectorTools::point_value(forward_dof_handler, forward_solution, p); + const double eps_n = Utilities::generate_normal_random_number(0, noise_level.back()); + + measurement_values.push_back (z_n + eps_n); + } + + std::cout << std::endl; +} + + +// It will make our lives easier if we can always assume that detector +// locations are at cell centers, because then we can evaluate the +// solution there using a quadrature formula whose sole quadrature +// point lies at the center of a cell. That's of course not where the +// "real" detector locations are, but it does not introduce a large +// error to do this. +template +void InformationDensityMeshRefinement::bounce_measurement_points_to_cell_centers () +{ + detector_locations_on_mesh = detector_locations; + for (auto &p : detector_locations_on_mesh) + { + for (const auto &cell : triangulation.active_cell_iterators()) + if (cell->point_inside (p)) + { + p = cell->center(); + break; + } + } +} + + +// The following functions are all quite standard by what we have +// shown in step-4, step-6, and step-22 (to name just a few of the +// more typical programs): +template +void InformationDensityMeshRefinement::setup_system () +{ + std::cout << "Setting up the linear system for the inverse problem..." + << std::endl; + + dof_handler.distribute_dofs (fe); + DoFRenumbering::component_wise (dof_handler); + + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints(dof_handler, + hanging_node_constraints); + hanging_node_constraints.close(); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + const std::vector dofs_per_component = + DoFTools::count_dofs_per_fe_component(dof_handler); + BlockDynamicSparsityPattern c_sparsity(dofs_per_component,dofs_per_component); + DoFTools::make_sparsity_pattern (dof_handler, c_sparsity); + hanging_node_constraints.condense(c_sparsity); + sparsity_pattern.copy_from(c_sparsity); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dofs_per_component); + system_rhs.reinit (dofs_per_component); +} + + + +template +void InformationDensityMeshRefinement::assemble_system () +{ + std::cout << "Assembling the linear system for the inverse problem..." + << std::flush; + + QGauss quadrature_formula(3); + + FEValues fe_values (fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + Vector cell_rhs (dofs_per_cell); + + std::vector local_dof_indices (dofs_per_cell); + + FEValuesExtractors::Scalar c(0), lambda(1), f(2); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit (cell); + cell_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point=0; q_point grad_phi_i = fe_values[c].gradient (i,q_point); + const Tensor<1,dim> grad_psi_i = fe_values[lambda].gradient (i,q_point); + + const double phi_i = fe_values[c].value (i,q_point); + const double psi_i = fe_values[lambda].value (i,q_point); + const double chi_i = fe_values[f].value (i,q_point); + + for (unsigned int j=0; j grad_phi_j = fe_values[c].gradient (j,q_point); + const Tensor<1,dim> grad_psi_j = fe_values[lambda].gradient (j,q_point); + + const double phi_j = fe_values[c].value (j,q_point); + const double psi_j= fe_values[lambda].value (j,q_point); + const double chi_j = fe_values[f].value (j,q_point); + + cell_matrix(i,j) += + ((grad_phi_i * grad_phi_j + + + phi_i * (velocity * grad_phi_j) + - + phi_i * chi_j + + + grad_psi_i * grad_psi_j + - + psi_i * (velocity * grad_psi_j) + - + chi_i * psi_j + + + regularization_parameter * chi_i * chi_j + ) * + fe_values.JxW (q_point)); + + for (unsigned int n=0; n< detector_locations_on_mesh.size(); ++n) + if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12) + { + cell_matrix(i,j) += psi_i * phi_j / noise_level[n] / noise_level[n]; + } + } + + for (unsigned int n=0; n< detector_locations_on_mesh.size(); ++n) + if (fe_values.quadrature_point(q_point).distance (detector_locations_on_mesh[n]) < 1e-12) + cell_rhs(i) += psi_i * measurement_values[n] / noise_level[n] / noise_level[n]; + } + + cell->get_dof_indices (local_dof_indices); + for (unsigned int i=0; i boundary_values; + std::vector component_mask (3); + component_mask[0] = component_mask[1] = true; + component_mask[2] = false; + VectorTools::interpolate_boundary_values (dof_handler, + 0, + Functions::ZeroFunction(3), + boundary_values, + component_mask); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + solution, + system_rhs); + + std::cout << std::endl; +} + + + +template +void InformationDensityMeshRefinement::solve () +{ + std::cout << "Solving the linear system for the inverse problem..." + << std::flush; + + SparseDirectUMFPACK A_direct; + solution = system_rhs; + A_direct.solve(system_matrix, solution); + + hanging_node_constraints.distribute (solution); + + std::cout << std::endl; +} + + + +// This is really the only interesting function of this program. It +// computes the functions $h_K = A^{-1} s_K$ for each source function +// (corresponding to each cell of the mesh). To do so, it first +// computes the forward matrix $A$ and uses the SparseDirectUMFPACK +// class to build an LU decomposition for this matrix. Then it loops +// over all cells $K$ and computes the corresponding $h_K$ by applying +// the LU decomposition to a right hand side vector for each $s_K$. +// +// The actual information content is then computed by evaluating these +// functions $h_K$ at measurement locations. +template +void InformationDensityMeshRefinement::compute_information_content () +{ + std::cout << "Computing the information content..." + << std::flush; + + information_content.reinit (triangulation.n_active_cells()); + + const FE_Q information_fe (fe.base_element(0).degree); + DoFHandler information_dof_handler (triangulation); + information_dof_handler.distribute_dofs (information_fe); + + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(information_dof_handler, constraints); + constraints.close(); + + SparsityPattern sparsity (information_dof_handler.n_dofs(), + information_dof_handler.max_couplings_between_dofs()); + DoFTools::make_sparsity_pattern (information_dof_handler, sparsity); + constraints.condense (sparsity); + sparsity.compress (); + + SparseMatrix system_matrix (sparsity); + + QGauss quadrature_formula(3); + + const unsigned int dofs_per_cell = information_fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + + // First build the forward operator + { + FEValues fe_values (information_fe, quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + for (const auto &cell : information_dof_handler.active_cell_iterators()) + { + fe_values.reinit (cell); + cell_matrix = 0; + + for (unsigned int q_point=0; q_pointdistribute_local_to_global (cell_matrix, + system_matrix); + } + + constraints.condense (system_matrix); + + std::map boundary_values; + VectorTools::interpolate_boundary_values (information_dof_handler, + 0, + Functions::ZeroFunction(), + boundary_values); + Vector tmp (information_dof_handler.n_dofs()); + MatrixTools::apply_boundary_values (boundary_values, + system_matrix, + tmp, + tmp); + } + + // Then factorize + SparseDirectUMFPACK A_inverse; + A_inverse.factorize(system_matrix); + + // Now compute the solutions corresponding to the possible + // sources. Each source is active on exactly one cell. + // + // As mentioned in the paper, this is a trivially parallel job, so + // we send the computations for each of these cells onto a separate + // task and let the OS schedule them onto individual processor + // cores. + Threads::TaskGroup tasks; + for (unsigned int K=0; K rhs (information_dof_handler.n_dofs()); + Vector cell_rhs (dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = information_dof_handler.begin_active(); + std::advance (cell, K); + + FEValues fe_values (information_fe, quadrature_formula, + update_values | + update_quadrature_points | update_JxW_values); + + fe_values.reinit (cell); + cell_rhs = 0; + + for (unsigned int q_point=0; q_pointdistribute_local_to_global (cell_rhs, + rhs); + + constraints.condense (rhs); + + A_inverse.solve(rhs); + + constraints.distribute (rhs); + + // Having computed the forward solutions + // corresponding to this source term, evaluate its + // contribution to the information content on all + // cells of the mesh by taking into account the + // detector locations. We add these into global + // objects, so we have to guard access to the + // global object: + static std::mutex m; + std::lock_guard g(m); + + + information_content(K) = regularization_parameter * cell->measure() * cell->measure(); + std::vector local_h_K_values (n_q_points); + for (const auto &cell : information_dof_handler.active_cell_iterators()) + { + fe_values.reinit (cell); + fe_values.get_function_values (rhs, local_h_K_values); + + for (unsigned int q_point=0; q_point +void InformationDensityMeshRefinement::output_results (const unsigned int cycle) const +{ + std::cout << "Outputting solutions..." << std::flush; + + DataOut data_out; + + std::vector names; + names.push_back ("forward_solution"); + names.push_back ("adjoint_solution"); + names.push_back ("recovered_parameter"); + + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, names); + data_out.add_data_vector (information_content, "information_content"); + + Vector information_density (triangulation.n_active_cells()); + for (const auto &cell : triangulation.active_cell_iterators()) + information_density(cell->active_cell_index()) + = std::sqrt(information_content(cell->active_cell_index())) / cell->measure(); + data_out.add_data_vector (information_density, "information_density"); + + data_out.build_patches (); + + std::string filename = "solution-"; + filename += ('0'+cycle); + filename += ".vtu"; + + std::ofstream output (filename.c_str()); + data_out.write_vtu (output); + + + // Now output the individual blocks of the matrix into files. + auto write_block = [&](const unsigned int block_i, + const unsigned int block_j, + const std::string &filename) + { + std::ofstream o(filename); + system_matrix.block(block_i,block_j).print (o); + }; + write_block(0,0, "matrix-" + std::to_string(cycle) + "-A.txt"); + write_block(0,2, "matrix-" + std::to_string(cycle) + "-B.txt"); + write_block(1,0, "matrix-" + std::to_string(cycle) + "-C.txt"); + write_block(2,2, "matrix-" + std::to_string(cycle) + "-M.txt"); + + std::cout << std::endl; +} + + + +// The following is then a function that refines the mesh based on the +// refinement criteria described in the paper. Which criterion to use +// is determined by which value the `refinement_criterion` variable +// is set to. +template +void InformationDensityMeshRefinement::refine_grid () +{ + std::cout << "Refining the mesh..." << std::endl; + + enum RefinementCriterion + { + global, + information_content, + indicator, + smoothness + }; + const RefinementCriterion refinement_criterion = information_content; + + switch (refinement_criterion) + { + case global: + { + triangulation.refine_global(); + break; + } + + case information_content: + { + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + this->information_content, + 0.2, 0.05); + triangulation.execute_coarsening_and_refinement (); + break; + } + + case indicator: + { + Vector refinement_indicators (triangulation.n_active_cells()); + + QGauss quadrature(3); + FEValues fe_values (fe, quadrature, update_values | update_JxW_values); + + FEValuesExtractors::Scalar lambda(1), f(2); + + std::vector lambda_values (quadrature.size()); + std::vector f_values (quadrature.size()); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit (cell); + fe_values[lambda].get_function_values (solution, lambda_values); + fe_values[f].get_function_values (solution, f_values); + + for (unsigned int q=0; qactive_cell_index()) + += (std::fabs (regularization_parameter * f_values[q] + - + lambda_values[q]) + * fe_values.JxW(q)); + } + + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + refinement_indicators, + 0.2, 0.05); + triangulation.execute_coarsening_and_refinement (); + break; + } + + + case smoothness: + { + Vector refinement_indicators (triangulation.n_active_cells()); + + DerivativeApproximation::approximate_gradient(dof_handler, + solution, + refinement_indicators, + /*component=*/2); + // and scale it to obtain an error indicator. + for (const auto &cell : triangulation.active_cell_iterators()) + refinement_indicators[cell->active_cell_index()] *= + std::pow(cell->diameter(), 1 + 1.0 * dim / 2); + + + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + refinement_indicators, + 0.2, 0.05); + triangulation.execute_coarsening_and_refinement (); + break; + } + + default: + Assert (false, ExcInternalError()); + } + + bounce_measurement_points_to_cell_centers (); + + + std::cout << std::endl; +} + + + + +template +void InformationDensityMeshRefinement::run () +{ + std::cout << "Solving problem in " << dim << " space dimensions." << std::endl; + + compute_synthetic_measurements (); + bounce_measurement_points_to_cell_centers (); + + for (unsigned int cycle=0; cycle<7; ++cycle) + { + std::cout << "---------- Cycle " << cycle << " ------------" << std::endl; + + setup_system (); + assemble_system (); + solve (); + compute_information_content (); + output_results (cycle); + refine_grid (); + } +} + + + +int main () +{ + try + { + deallog.depth_console (0); + + InformationDensityMeshRefinement<2> information_density_mesh_refinement; + information_density_mesh_refinement.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +}