coupled_laplace_problem.h

Go to the documentation of this file.

) const
 * {
 *   double return_value = 0.0;
 *   for (unsigned int i = 0; i < dim; ++i)
 *     return_value += 4.0 * std::pow(p(i), 4.0);
 * 
 *   return return_value;
 * }
 * 
 * 
 * template <int dim>
 * double
 * BoundaryValues<dim>::value(const Point<dim> &p,
 *                            const unsigned int /*component*/) const
 * {
 *   return p.square();
 * }
 * 
 * 
 * 
 * template <int dim>
 * CoupledLaplaceProblem<dim>::CoupledLaplaceProblem()
 *   : fe(1)
 *   , dof_handler(triangulation)
 *   , interface_boundary_id(1)
 *   , adapter(parameters, interface_boundary_id)
 * {}
 * 
 * 
 * template <int dim>
 * void
 * CoupledLaplaceProblem<dim>::make_grid()
 * {
 *   GridGenerator::hyper_cube(triangulation, -1, 1);
 *   triangulation.refine_global(4);
 * 
 *   for (const auto &cell : triangulation.active_cell_iterators())
 *     for (const auto face : cell->face_iterators())
 *       {
 * @endcode
 * 
 * We choose the boundary in positive x direction for the
 * interface coupling.
 * 
 * @code
 *         if (face->at_boundary() && (face->center()[0] == 1))
 *           face->set_boundary_id(interface_boundary_id);
 *       }
 * 
 *   std::cout << "   Number of active cells: " << triangulation.n_active_cells()
 *             << std::endl
 *             << "   Total number of cells: " << triangulation.n_cells()
 *             << std::endl;
 * }
 * 
 * 
 * template <int dim>
 * void
 * CoupledLaplaceProblem<dim>::setup_system()
 * {
 *   dof_handler.distribute_dofs(fe);
 * 
 *   std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
 *             << std::endl;
 * 
 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *   sparsity_pattern.copy_from(dsp);
 * 
 *   system_matrix.reinit(sparsity_pattern);
 * 
 *   solution.reinit(dof_handler.n_dofs());
 *   old_solution.reinit(dof_handler.n_dofs());
 *   system_rhs.reinit(dof_handler.n_dofs());
 * }
 * 
 * 
 * 
 * template <int dim>
 * void
 * CoupledLaplaceProblem<dim>::assemble_system()
 * {
 * @endcode
 * 
 * Reset global structures
 * 
 * @code
 *   system_rhs    = 0;
 *   system_matrix = 0;
 * @endcode
 * 
 * Update old solution values
 * 
 * @code
 *   old_solution = solution;
 * 
 *   QGauss<dim> quadrature_formula(fe.degree + 1);
 * 
 *   RightHandSide<dim> right_hand_side;
 * 
 *   FEValues<dim> fe_values(fe,
 *                           quadrature_formula,
 *                           update_values | update_gradients |
 *                             update_quadrature_points | update_JxW_values);
 * 
 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 * 
 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *   Vector<double>     cell_rhs(dofs_per_cell);
 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 * @endcode
 * 
 * The solution values from previous time steps are stored for each quadrature
 * point
 * 
 * @code
 *   std::vector<double> local_values_old_solution(fe_values.n_quadrature_points);
 * 
 *   for (const auto &cell : dof_handler.active_cell_iterators())
 *     {
 *       fe_values.reinit(cell);
 *       cell_matrix = 0;
 *       cell_rhs    = 0;
 * @endcode
 * 
 * Get the local values from the `fe_values' object
 * 
 * @code
 *       fe_values.get_function_values(old_solution, local_values_old_solution);
 * 
 * @endcode
 * 
 * The system matrix contains additionally a mass matrix due to the time
 * discretization. The RHS has contributions from the old solution values.
 * 
 * @code
 *       for (const unsigned int q_index : fe_values.quadrature_point_indices())
 *         for (const unsigned int i : fe_values.dof_indices())
 *           {
 *             for (const unsigned int j : fe_values.dof_indices())
 *               cell_matrix(i, j) +=
 *                 ((fe_values.shape_value(i, q_index) *  // phi_i(x_q)
 *                   fe_values.shape_value(j, q_index)) + // phi_j(x_q)
 *                  (delta_t *                            // delta t
 *                   fe_values.shape_grad(i, q_index) *   // grad phi_i(x_q)
 *                   fe_values.shape_grad(j, q_index))) * // grad phi_j(x_q)
 *                 fe_values.JxW(q_index);                // dx
 * 
 *             const auto  x_q         = fe_values.quadrature_point(q_index);
 *             const auto &local_value = local_values_old_solution[q_index];
 *             cell_rhs(i) += ((delta_t *                           // delta t
 *                              fe_values.shape_value(i, q_index) * // phi_i(x_q)
 *                              right_hand_side.value(x_q)) +       // f(x_q)
 *                             fe_values.shape_value(i, q_index) *
 *                               local_value) *       // phi_i(x_q)*val
 *                            fe_values.JxW(q_index); // dx
 *           }
 * 
 * @endcode
 * 
 * Copy local to global
 * 
 * @code
 *       cell->get_dof_indices(local_dof_indices);
 *       for (const unsigned int i : fe_values.dof_indices())
 *         {
 *           for (const unsigned int j : fe_values.dof_indices())
 *             system_matrix.add(local_dof_indices[i],
 *                               local_dof_indices[j],
 *                               cell_matrix(i, j));
 * 
 *           system_rhs(local_dof_indices[i]) += cell_rhs(i);
 *         }
 *     }
 *   {
 * @endcode
 * 
 * At first, we apply the Dirichlet boundary condition from @ref step_4 "step-4", as
 * usual.
 * 
 * @code
 *     std::map<types::global_dof_index, double> boundary_values;
 *     VectorTools::interpolate_boundary_values(dof_handler,
 *                                              0,
 *                                              BoundaryValues<dim>(),
 *                                              boundary_values);
 *     MatrixTools::apply_boundary_values(boundary_values,
 *                                        system_matrix,
 *                                        solution,
 *                                        system_rhs);
 *   }
 *   {
 * @endcode
 * 
 * Afterwards, we apply the coupling boundary condition. The `boundary_data`
 * has already been filled by preCICE.
 * 
 * @code
 *     MatrixTools::apply_boundary_values(boundary_data,
 *                                        system_matrix,
 *                                        solution,
 *                                        system_rhs);
 *   }
 * }
 * 
 * 
 * 
 * template <int dim>
 * void
 * CoupledLaplaceProblem<dim>::solve()
 * {
 *   SolverControl            solver_control(1000, 1e-12);
 *   SolverCG<Vector<double>> solver(solver_control);
 *   solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
 * 
 *   std::cout << "   " << solver_control.last_step()
 *             << " CG iterations needed to obtain convergence." << std::endl;
 * }
 * 
 * 
 * 
 * template <int dim>
 * void
 * CoupledLaplaceProblem<dim>::output_results() const
 * {
 *   DataOut<dim> data_out;
 * 
 *   data_out.attach_dof_handler(dof_handler);
 *   data_out.add_data_vector(solution, "solution");
 * 
 *   data_out.build_patches(mapping);
 * 
 *   std::ofstream output("solution-" + std::to_string(time_step) + ".vtk");
 *   data_out.write_vtk(output);
 * }
 * 
 * 
 * 
 * template <int dim>
 * void
 * CoupledLaplaceProblem<dim>::run()
 * {
 *   std::cout << "Solving problem in " << dim << " space dimensions."
 *             << std::endl;
 * 
 *   make_grid();
 *   setup_system();
 * 
 * @endcode
 * 
 * After we set up the system, we initialize preCICE using the functionalities
 * of the Adapter. preCICE returns the maximum admissible time-step size,
 * which needs to be compared to our desired solver time-step size.
 * 
 * @code
 *   precice_delta_t = adapter.initialize(dof_handler, boundary_data, mapping);
 *   delta_t         = std::min(precice_delta_t, solver_delta_t);
 * 
 * @endcode
 * 
 * preCICE steers the coupled simulation: `isCouplingOngoing` is
 * used to synchronize the end of the simulation with the coupling partner
 * 
 * @code
 *   while (adapter.precice.isCouplingOngoing())
 *     {
 * @endcode
 * 
 * The time step number is solely used to generate unique output files
 * 
 * @code
 *       ++time_step;
 * @endcode
 * 
 * In the time loop, we assemble the coupled system and solve it as
 * usual.
 * 
 * @code
 *       assemble_system();
 *       solve();
 * 
 * @endcode
 * 
 * After we solved the system, we advance the coupling to the next time
 * level. In a bi-directional coupled simulation, we would pass our
 * calculated data to and obtain new data from preCICE. Here, we simply
 * obtain new data from preCICE, so from the other participant. As before,
 * we obtain a maximum time-step size and compare it against the desired
 * solver time-step size.
 * 
 * @code
 *       precice_delta_t = adapter.advance(boundary_data, delta_t);
 *       delta_t         = std::min(precice_delta_t, solver_delta_t);
 * 
 * @endcode
 * 
 * Write an output file if the time step is completed. In case of an
 * implicit coupling, where individual time steps are computed more than
 * once, the function `isTimeWindowCompleted` prevents unnecessary result
 * writing. For this simple tutorial configuration (explicit coupling),
 * the function returns always `true`.
 * 
 * @code
 *       if (adapter.precice.isTimeWindowComplete())
 *         output_results();
 *     }
 * }
 * 
 * 
 * 
 * int
 * main()
 * {
 *   CoupledLaplaceProblem<2> laplace_problem;
 *   laplace_problem.run();
 * 
 *   return 0;
 * }
 * @endcode
 
 
<a name="ann-fancy_boundary_condition.cc"></a>
<h1>Annotated version of fancy_boundary_condition.cc</h1>
 * 
 * 
 * 
 * This program does not use any deal.II functionality and depends only on
 * preCICE and the standard libraries.
 * 
 * @code
 * #include <precice/SolverInterface.hpp>
 * 
 * #include <iostream>
 * #include <sstream>
 * 
 * @endcode
 * 
 * The program computes a time-varying parabolic boundary condition, which is
 * passed to preCICE and serves as Dirichlet boundary condition for the other
 * coupling participant.
 * 
 
 * 
 * Function to generate boundary values in each time step
 * 
 * @code
 * void
 * define_boundary_values(std::vector<double> &boundary_data,
 *                        const double         time,
 *                        const double         end_time)
 * {
 * @endcode
 * 
 * Scale the current time value
 * 
 * @code
 *   const double relative_time = time / end_time;
 * @endcode
 * 
 * Define the amplitude. Values run from -0.5 to 0.5
 * 
 * @code
 *   const double amplitude = (relative_time - 0.5);
 * @endcode
 * 
 * Specify the actual data we want to pass to the other participant. Here, we
 * choose a parabola with boundary values 2 in order to enforce continuity
 * to adjacent boundaries.
 * 
 * @code
 *   const double n_elements = boundary_data.size();
 *   const double right_zero = boundary_data.size() - 1;
 *   const double left_zero  = 0;
 *   const double offset     = 2;
 *   for (uint i = 0; i < n_elements; ++i)
 *     boundary_data[i] =
 *       -amplitude * ((i - left_zero) * (i - right_zero)) + offset;
 * }
 * 
 * 
 * int
 * main()
 * {
 *   std::cout << "Boundary participant: starting... \n";
 * 
 * @endcode
 * 
 * Configuration
 * 
 * @code
 *   const std::string configFileName("precice-config.xml");
 *   const std::string solverName("boundary-participant");
 *   const std::string meshName("boundary-mesh");
 *   const std::string dataWriteName("boundary-data");
 * 
 * @endcode
 * 
 * Adjust to MPI rank and size for parallel computation
 * 
 * @code
 *   const int commRank = 0;
 *   const int commSize = 1;
 * 
 *   precice::SolverInterface precice(solverName,
 *                                    configFileName,
 *                                    commRank,
 *                                    commSize);
 * 
 *   const int meshID           = precice.getMeshID(meshName);
 *   const int dimensions       = precice.getDimensions();
 *   const int numberOfVertices = 6;
 * 
 *   const int dataID = precice.getDataID(dataWriteName, meshID);
 * 
 * @endcode
 * 
 * Set up data structures
 * 
 * @code
 *   std::vector<double> writeData(numberOfVertices);
 *   std::vector<int>    vertexIDs(numberOfVertices);
 *   std::vector<double> vertices(numberOfVertices * dimensions);
 * 
 * @endcode
 * 
 * Define a boundary mesh
 * 
 * @code
 *   std::cout << "Boundary participant: defining boundary mesh \n";
 *   const double length = 2;
 *   const double xCoord = 1;
 *   const double deltaY = length / (numberOfVertices - 1);
 *   for (int i = 0; i < numberOfVertices; ++i)
 *     for (int j = 0; j < dimensions; ++j)
 *       {
 *         const unsigned int index = dimensions * i + j;
 * @endcode
 * 
 * The x-coordinate is always 1, i.e., the boundary is parallel to the
 * y-axis. The y-coordinate is descending from 1 to -1.
 * 
 * @code
 *         if (j == 0)
 *           vertices[index] = xCoord;
 *         else
 *           vertices[index] = 1 - deltaY * i;
 *       }
 * 
 * @endcode
 * 
 * Pass the vertices to preCICE
 * 
 * @code
 *   precice.setMeshVertices(meshID,
 *                           numberOfVertices,
 *                           vertices.data(),
 *                           vertexIDs.data());
 * 
 * @endcode
 * 
 * initialize the Solverinterface
 * 
 * @code
 *   double dt = precice.initialize();
 * 
 * @endcode
 * 
 * Start time loop
 * 
 * @code
 *   const double end_time = 1;
 *   double       time     = 0;
 *   while (precice.isCouplingOngoing())
 *     {
 * @endcode
 * 
 * Generate new boundary data
 * 
 * @code
 *       define_boundary_values(writeData, time, end_time);
 * 
 *       {
 *         std::cout << "Boundary participant: writing coupling data \n";
 *         precice.writeBlockScalarData(dataID,
 *                                      numberOfVertices,
 *                                      vertexIDs.data(),
 *                                      writeData.data());
 *       }
 * 
 *       dt = precice.advance(dt);
 *       std::cout << "Boundary participant: advancing in time\n";
 * 
 *       time += dt;
 *     }
 * 
 *   std::cout << "Boundary participant: closing...\n";
 * 
 *   return 0;
 * }
 * @endcode
 
 
*/