198 *
void refine_grid ();
248 *
for (
unsigned int i=0; i<50; ++i)
260 *
for (
unsigned int i=0; i<50; ++i)
303 *
std::ofstream output(
"detector_locations.vtu");
309 *
While we're generating output, also output the source location. Do this
310 * by outputting many (1000) points that indicate the perimeter of the source
314 * Particles::ParticleHandler<dim> particle_handler(triangulation,
315 * StaticMappingQ1<dim>::mapping);
317 * const unsigned int n_points = 1000;
318 * for (unsigned int i=0; i<n_points; ++i)
320 * Point<dim> loc = source_location;
321 * loc[0] += source_radius * std::cos(2*numbers::PI*i/n_points);
322 * loc[1] += source_radius * std::sin(2*numbers::PI*i/n_points);
324 * Particles::Particle<dim> new_particle;
325 * new_particle.set_location(loc);
326 * particle_handler.insert_particle(new_particle,
327 * triangulation.begin_active());
330 * Particles::DataOut<dim> particle_out;
331 * particle_out.build_patches(particle_handler);
332 * std::ofstream output("source_locations.vtu");
333 * particle_out.write_vtu(output);
341 * The following function solves a forward problem on a twice
342 * refined mesh to compute "synthetic data". Refining the mesh
343 * beyond the mesh used for the inverse problem avoids an
348 * void InformationDensityMeshRefinement<dim>::compute_synthetic_measurements ()
350 * std::cout << "Computing synthetic data by solving the forward problem..."
355 * Create a triangulation and DoFHandler that corresponds to a
356 * twice-refined mesh so that we obtain the synthetic data with
357 * higher accuracy than we do on the regular mesh used for all other
361 * Triangulation<dim> forward_triangulation;
362 * forward_triangulation.copy_triangulation (triangulation);
363 * forward_triangulation.refine_global (2);
365 * const FE_Q<dim> forward_fe (fe.base_element(0).degree);
366 * DoFHandler<dim> forward_dof_handler (forward_triangulation);
367 * forward_dof_handler.distribute_dofs (forward_fe);
369 * AffineConstraints<double> constraints;
370 * DoFTools::make_hanging_node_constraints(forward_dof_handler, constraints);
371 * constraints.close();
373 * SparsityPattern sparsity (forward_dof_handler.n_dofs(),
374 * forward_dof_handler.max_couplings_between_dofs());
375 * DoFTools::make_sparsity_pattern (forward_dof_handler, sparsity);
376 * constraints.condense (sparsity);
377 * sparsity.compress ();
379 * SparseMatrix<double> system_matrix (sparsity);
380 * Vector<double> system_rhs (forward_dof_handler.n_dofs());
382 * QGauss<dim> quadrature_formula(3);
383 * FEValues<dim> fe_values (forward_fe, quadrature_formula,
384 * update_values | update_gradients |
385 * update_quadrature_points | update_JxW_values);
387 * const unsigned int dofs_per_cell = forward_fe.dofs_per_cell;
388 * const unsigned int n_q_points = quadrature_formula.size();
392 * First assemble the system matrix and right hand side for the forward
397 * FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
398 * Vector<double> cell_rhs (dofs_per_cell);
399 * std::vector<unsigned int> local_dof_indices (dofs_per_cell);
401 * for (const auto &cell : forward_dof_handler.active_cell_iterators())
403 * fe_values.reinit (cell);
407 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
408 * for (unsigned int i=0; i<dofs_per_cell; ++i)
409 * for (unsigned int j=0; j<dofs_per_cell; ++j)
410 * cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
411 * fe_values.shape_grad(j,q_point)
413 * fe_values.shape_value(i,q_point) *
414 * (velocity * fe_values.shape_grad(j,q_point))
417 * fe_values.JxW(q_point);
418 * for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
419 * if (fe_values.quadrature_point(q_point).distance (source_location)
421 * for (unsigned int i=0; i<dofs_per_cell; ++i)
424 * fe_values.shape_value (i, q_point) *
425 * fe_values.JxW(q_point);
427 * cell->get_dof_indices (local_dof_indices);
428 * constraints.distribute_local_to_global (cell_matrix,
435 * std::map<unsigned int, double> boundary_values;
436 * VectorTools::interpolate_boundary_values (forward_dof_handler,
438 * Functions::ZeroFunction<dim>(),
440 * Vector<double> tmp (forward_dof_handler.n_dofs());
441 * MatrixTools::apply_boundary_values (boundary_values,
449 * Solve the forward problem and output it into its own VTU file :
452 * SparseDirectUMFPACK A_inverse;
453 * Vector<double> forward_solution (forward_dof_handler.n_dofs());
454 * forward_solution = system_rhs;
455 * A_inverse.solve(system_matrix, forward_solution);
457 * const double max_forward_solution = forward_solution.linfty_norm();
460 * DataOut<dim> data_out;
461 * data_out.attach_dof_handler (forward_dof_handler);
462 * data_out.add_data_vector (forward_solution, "c");
463 * data_out.build_patches (4);
465 * std::ofstream out ("forward-solution.vtu");
466 * data_out.write_vtu (out);
471 * Now evaluate the forward solution at the measurement points:
474 * for (const auto &p : detector_locations)
478 * same 10% noise level for all points
481 * noise_level.push_back (0.1 * max_forward_solution);
483 * const double z_n = VectorTools::point_value(forward_dof_handler, forward_solution, p);
484 * const double eps_n = Utilities::generate_normal_random_number(0, noise_level.back());
486 * measurement_values.push_back (z_n + eps_n);
489 * std::cout << std::endl;
495 * It will make our lives easier if we can always assume that detector
496 * locations are at cell centers, because then we can evaluate the
497 * solution there using a quadrature formula whose sole quadrature
510 *
if (cell->point_inside (p))
512 *
p = cell->center();
529 *
std::
cout <<
"Setting up the linear system for the inverse problem..."
532 *
dof_handler.distribute_dofs (fe);
540 *
std::cout <<
" Number of active cells: "
543 *
std::cout <<
" Number of degrees of freedom: "
544 *
<< dof_handler.n_dofs()
547 *
const std::vector<types::global_dof_index> dofs_per_component =
554 *
system_matrix.reinit (sparsity_pattern);
556 *
solution.reinit (dofs_per_component);
565 *
std::cout <<
"Assembling the linear system for the inverse problem..."
574 *
const unsigned int dofs_per_cell = fe.dofs_per_cell;
580 *
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
584 *
for (
const auto &cell : dof_handler.active_cell_iterators())
586 *
fe_values.
reinit (cell);
591 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
600 *
for (
unsigned int j=0;
j<dofs_per_cell; ++
j)
638 *
cell->get_dof_indices (local_dof_indices);
639 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
641 *
for (
unsigned int j=0;
j<dofs_per_cell; ++
j)
642 *
system_matrix.add (local_dof_indices[i],
643 *
local_dof_indices[
j],
654 *
std::vector<bool> component_mask (3);
655 *
component_mask[0] = component_mask[1] =
true;
656 *
component_mask[2] =
false;
667 *
std::cout << std::endl;
675 *
std::cout <<
"Solving the linear system for the inverse problem..."
680 *
A_direct.solve(system_matrix, solution);
684 *
std::cout << std::endl;
708 *
std::cout <<
"Computing the information content..."
719 *
constraints.close();
724 *
constraints.condense (sparsity);
725 *
sparsity.compress ();
731 *
const unsigned int dofs_per_cell =
information_fe.dofs_per_cell;
745 *
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
749 *
fe_values.
reinit (cell);
753 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
754 *
for (
unsigned int j=0;
j<dofs_per_cell; ++
j)
758 *
fe_values.shape_value(i,
q_point) *
762 *
cell->distribute_local_to_global (cell_matrix,
766 *
constraints.condense (system_matrix);
808 *
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
812 *
std::advance (cell, K);
818 *
fe_values.reinit (cell);
822 *
for (
unsigned int i=0; i<dofs_per_cell; ++i)
826 *
cell->distribute_local_to_global (
cell_rhs,
829 *
constraints.condense (
rhs);
833 *
constraints.distribute (
rhs);
846 *
static std::mutex m;
847 *
std::lock_guard<std::mutex>
g(m);
854 *
fe_values.
reinit (cell);
874 *
std::cout << std::endl;
890 *
std::cout <<
"Outputting solutions..." << std::flush;
894 *
std::vector<std::string> names;
895 *
names.push_back (
"forward_solution");
896 *
names.push_back (
"adjoint_solution");
897 *
names.push_back (
"recovered_parameter");
899 *
data_out.attach_dof_handler (dof_handler);
900 *
data_out.add_data_vector (solution, names);
904 *
for (
const auto &cell :
triangulation.active_cell_iterators())
909 *
data_out.build_patches ();
911 *
std::string
filename =
"solution-";
915 *
std::ofstream output (
filename.c_str());
916 *
data_out.write_vtu (output);
931 *
write_block(0,0,
"matrix-" + std::to_string(cycle) +
"-A.txt");
932 *
write_block(0,2,
"matrix-" + std::to_string(cycle) +
"-B.txt");
933 *
write_block(1,0,
"matrix-" + std::to_string(cycle) +
"-C.txt");
934 *
write_block(2,2,
"matrix-" + std::to_string(cycle) +
"-M.txt");
936 *
std::cout << std::endl;
952 *
std::cout <<
"Refining the mesh..." << std::endl;
990 *
std::vector<double>
f_values (quadrature.size());
992 *
for (
const auto &cell : dof_handler.active_cell_iterators())
994 *
fe_values.
reinit (cell);
996 *
fe_values[f].get_function_values (solution,
f_values);
998 *
for (
unsigned int q=0;
q<quadrature.size(); ++
q)
1003 *
* fe_values.JxW(
q));
1027 *
for (
const auto &cell :
triangulation.active_cell_iterators())
1046 *
std::cout << std::endl;
1052 *
template <
int dim>
1055 *
std::cout <<
"Solving problem in " << dim <<
" space dimensions." << std::endl;
1060 *
for (
unsigned int cycle=0; cycle<7; ++cycle)
1062 *
std::cout <<
"---------- Cycle " << cycle <<
" ------------" << std::endl;
1084 *
catch (std::exception &exc)
1086 *
std::cerr << std::endl << std::endl
1087 *
<<
"----------------------------------------------------"
1089 *
std::cerr <<
"Exception on processing: " << std::endl
1090 *
<< exc.what() << std::endl
1091 *
<<
"Aborting!" << std::endl
1092 *
<<
"----------------------------------------------------"
1099 *
std::cerr << std::endl << std::endl
1100 *
<<
"----------------------------------------------------"
1102 *
std::cerr <<
"Unknown exception!" << std::endl
1103 *
<<
"Aborting!" << std::endl
1104 *
<<
"----------------------------------------------------"
unsigned int depth_console(const unsigned int n)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
typename ActiveSelector::active_cell_iterator active_cell_iterator
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
Task< RT > new_task(const std::function< RT()> &function)
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void approximate_gradient(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InputVector &solution, Vector< float > &derivative_norm, const unsigned int component=0)
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
constexpr types::blas_int one
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)