Reference documentation for deal.II version 9.5.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
coupled_laplace_problem.h
Go to the documentation of this file.
1) const
651 *   {
652 *   double return_value = 0.0;
653 *   for (unsigned int i = 0; i < dim; ++i)
654 *   return_value += 4.0 * std::pow(p(i), 4.0);
655 *  
656 *   return return_value;
657 *   }
658 *  
659 *  
660 *   template <int dim>
661 *   double
662 *   BoundaryValues<dim>::value(const Point<dim> &p,
663 *   const unsigned int /*component*/) const
664 *   {
665 *   return p.square();
666 *   }
667 *  
668 *  
669 *  
670 *   template <int dim>
671 *   CoupledLaplaceProblem<dim>::CoupledLaplaceProblem()
672 *   : fe(1)
673 *   , dof_handler(triangulation)
674 *   , interface_boundary_id(1)
675 *   , adapter(parameters, interface_boundary_id)
676 *   {}
677 *  
678 *  
679 *   template <int dim>
680 *   void
681 *   CoupledLaplaceProblem<dim>::make_grid()
682 *   {
684 *   triangulation.refine_global(4);
685 *  
686 *   for (const auto &cell : triangulation.active_cell_iterators())
687 *   for (const auto &face : cell->face_iterators())
688 *   {
689 * @endcode
690 *
691 * We choose the boundary in positive x direction for the
692 * interface coupling.
693 *
694 * @code
695 *   if (face->at_boundary() && (face->center()[0] == 1))
696 *   face->set_boundary_id(interface_boundary_id);
697 *   }
698 *  
699 *   std::cout << " Number of active cells: " << triangulation.n_active_cells()
700 *   << std::endl
701 *   << " Total number of cells: " << triangulation.n_cells()
702 *   << std::endl;
703 *   }
704 *  
705 *  
706 *   template <int dim>
707 *   void
708 *   CoupledLaplaceProblem<dim>::setup_system()
709 *   {
710 *   dof_handler.distribute_dofs(fe);
711 *  
712 *   std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
713 *   << std::endl;
714 *  
715 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
716 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
717 *   sparsity_pattern.copy_from(dsp);
718 *  
719 *   system_matrix.reinit(sparsity_pattern);
720 *  
721 *   solution.reinit(dof_handler.n_dofs());
722 *   old_solution.reinit(dof_handler.n_dofs());
723 *   system_rhs.reinit(dof_handler.n_dofs());
724 *   }
725 *  
726 *  
727 *  
728 *   template <int dim>
729 *   void
730 *   CoupledLaplaceProblem<dim>::assemble_system()
731 *   {
732 * @endcode
733 *
734 * Reset global structures
735 *
736 * @code
737 *   system_rhs = 0;
738 *   system_matrix = 0;
739 * @endcode
740 *
741 * Update old solution values
742 *
743 * @code
744 *   old_solution = solution;
745 *  
746 *   QGauss<dim> quadrature_formula(fe.degree + 1);
747 *  
748 *   RightHandSide<dim> right_hand_side;
749 *  
750 *   FEValues<dim> fe_values(fe,
751 *   quadrature_formula,
754 *  
755 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
756 *  
757 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
758 *   Vector<double> cell_rhs(dofs_per_cell);
759 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
760 * @endcode
761 *
762 * The solution values from previous time steps are stored for each quadrature
763 * point
764 *
765 * @code
766 *   std::vector<double> local_values_old_solution(fe_values.n_quadrature_points);
767 *  
768 *   for (const auto &cell : dof_handler.active_cell_iterators())
769 *   {
770 *   fe_values.reinit(cell);
771 *   cell_matrix = 0;
772 *   cell_rhs = 0;
773 * @endcode
774 *
775 * Get the local values from the `fe_values' object
776 *
777 * @code
778 *   fe_values.get_function_values(old_solution, local_values_old_solution);
779 *  
780 * @endcode
781 *
782 * The system matrix contains additionally a mass matrix due to the time
783 * discretization. The RHS has contributions from the old solution values.
784 *
785 * @code
786 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
787 *   for (const unsigned int i : fe_values.dof_indices())
788 *   {
789 *   for (const unsigned int j : fe_values.dof_indices())
790 *   cell_matrix(i, j) +=
791 *   ((fe_values.shape_value(i, q_index) * // phi_i(x_q)
792 *   fe_values.shape_value(j, q_index)) + // phi_j(x_q)
793 *   (delta_t * // delta t
794 *   fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
795 *   fe_values.shape_grad(j, q_index))) * // grad phi_j(x_q)
796 *   fe_values.JxW(q_index); // dx
797 *  
798 *   const auto x_q = fe_values.quadrature_point(q_index);
799 *   const auto &local_value = local_values_old_solution[q_index];
800 *   cell_rhs(i) += ((delta_t * // delta t
801 *   fe_values.shape_value(i, q_index) * // phi_i(x_q)
802 *   right_hand_side.value(x_q)) + // f(x_q)
803 *   fe_values.shape_value(i, q_index) *
804 *   local_value) * // phi_i(x_q)*val
805 *   fe_values.JxW(q_index); // dx
806 *   }
807 *  
808 * @endcode
809 *
810 * Copy local to global
811 *
812 * @code
813 *   cell->get_dof_indices(local_dof_indices);
814 *   for (const unsigned int i : fe_values.dof_indices())
815 *   {
816 *   for (const unsigned int j : fe_values.dof_indices())
817 *   system_matrix.add(local_dof_indices[i],
818 *   local_dof_indices[j],
819 *   cell_matrix(i, j));
820 *  
821 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
822 *   }
823 *   }
824 *   {
825 * @endcode
826 *
827 * At first, we apply the Dirichlet boundary condition from @ref step_4 "step-4", as
828 * usual.
829 *
830 * @code
831 *   std::map<types::global_dof_index, double> boundary_values;
832 *   VectorTools::interpolate_boundary_values(dof_handler,
833 *   0,
834 *   BoundaryValues<dim>(),
835 *   boundary_values);
836 *   MatrixTools::apply_boundary_values(boundary_values,
837 *   system_matrix,
838 *   solution,
839 *   system_rhs);
840 *   }
841 *   {
842 * @endcode
843 *
844 * Afterwards, we apply the coupling boundary condition. The `boundary_data`
845 * has already been filled by preCICE.
846 *
847 * @code
848 *   MatrixTools::apply_boundary_values(boundary_data,
849 *   system_matrix,
850 *   solution,
851 *   system_rhs);
852 *   }
853 *   }
854 *  
855 *  
856 *  
857 *   template <int dim>
858 *   void
859 *   CoupledLaplaceProblem<dim>::solve()
860 *   {
861 *   SolverControl solver_control(1000, 1e-12);
862 *   SolverCG<Vector<double>> solver(solver_control);
863 *   solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
864 *  
865 *   std::cout << " " << solver_control.last_step()
866 *   << " CG iterations needed to obtain convergence." << std::endl;
867 *   }
868 *  
869 *  
870 *  
871 *   template <int dim>
872 *   void
873 *   CoupledLaplaceProblem<dim>::output_results() const
874 *   {
875 *   DataOut<dim> data_out;
876 *  
877 *   data_out.attach_dof_handler(dof_handler);
878 *   data_out.add_data_vector(solution, "solution");
879 *  
880 *   data_out.build_patches(mapping);
881 *  
882 *   std::ofstream output("solution-" + std::to_string(time_step) + ".vtk");
883 *   data_out.write_vtk(output);
884 *   }
885 *  
886 *  
887 *  
888 *   template <int dim>
889 *   void
890 *   CoupledLaplaceProblem<dim>::run()
891 *   {
892 *   std::cout << "Solving problem in " << dim << " space dimensions."
893 *   << std::endl;
894 *  
895 *   make_grid();
896 *   setup_system();
897 *  
898 * @endcode
899 *
900 * After we set up the system, we initialize preCICE using the functionalities
901 * of the Adapter. preCICE returns the maximum admissible time-step size,
902 * which needs to be compared to our desired solver time-step size.
903 *
904 * @code
905 *   precice_delta_t = adapter.initialize(dof_handler, boundary_data, mapping);
906 *   delta_t = std::min(precice_delta_t, solver_delta_t);
907 *  
908 * @endcode
909 *
910 * preCICE steers the coupled simulation: `isCouplingOngoing` is
911 * used to synchronize the end of the simulation with the coupling partner
912 *
913 * @code
914 *   while (adapter.precice.isCouplingOngoing())
915 *   {
916 * @endcode
917 *
918 * The time step number is solely used to generate unique output files
919 *
920 * @code
921 *   ++time_step;
922 * @endcode
923 *
924 * In the time loop, we assemble the coupled system and solve it as
925 * usual.
926 *
927 * @code
928 *   assemble_system();
929 *   solve();
930 *  
931 * @endcode
932 *
933 * After we solved the system, we advance the coupling to the next time
934 * level. In a bi-directional coupled simulation, we would pass our
935 * calculated data to and obtain new data from preCICE. Here, we simply
936 * obtain new data from preCICE, so from the other participant. As before,
937 * we obtain a maximum time-step size and compare it against the desired
938 * solver time-step size.
939 *
940 * @code
941 *   precice_delta_t = adapter.advance(boundary_data, delta_t);
942 *   delta_t = std::min(precice_delta_t, solver_delta_t);
943 *  
944 * @endcode
945 *
946 * Write an output file if the time step is completed. In case of an
947 * implicit coupling, where individual time steps are computed more than
948 * once, the function `isTimeWindowCompleted` prevents unnecessary result
949 * writing. For this simple tutorial configuration (explicit coupling),
950 * the function returns always `true`.
951 *
952 * @code
953 *   if (adapter.precice.isTimeWindowComplete())
954 *   output_results();
955 *   }
956 *   }
957 *  
958 *  
959 *  
960 *   int
961 *   main()
962 *   {
963 *   CoupledLaplaceProblem<2> laplace_problem;
964 *   laplace_problem.run();
965 *  
966 *   return 0;
967 *   }
968 * @endcode
969
970
971<a name="ann-fancy_boundary_condition.cc"></a>
972<h1>Annotated version of fancy_boundary_condition.cc</h1>
973 *
974 *
975 *
976 * This program does not use any deal.II functionality and depends only on
977 * preCICE and the standard libraries.
978 *
979 * @code
980 *   #include <precice/SolverInterface.hpp>
981 *  
982 *   #include <iostream>
983 *   #include <sstream>
984 *  
985 * @endcode
986 *
987 * The program computes a time-varying parabolic boundary condition, which is
988 * passed to preCICE and serves as Dirichlet boundary condition for the other
989 * coupling participant.
990 *
991
992 *
993 * Function to generate boundary values in each time step
994 *
995 * @code
996 *   void
997 *   define_boundary_values(std::vector<double> &boundary_data,
998 *   const double time,
999 *   const double end_time)
1000 *   {
1001 * @endcode
1002 *
1003 * Scale the current time value
1004 *
1005 * @code
1006 *   const double relative_time = time / end_time;
1007 * @endcode
1008 *
1009 * Define the amplitude. Values run from -0.5 to 0.5
1010 *
1011 * @code
1012 *   const double amplitude = (relative_time - 0.5);
1013 * @endcode
1014 *
1015 * Specify the actual data we want to pass to the other participant. Here, we
1016 * choose a parabola with boundary values 2 in order to enforce continuity
1017 * to adjacent boundaries.
1018 *
1019 * @code
1020 *   const double n_elements = boundary_data.size();
1021 *   const double right_zero = boundary_data.size() - 1;
1022 *   const double left_zero = 0;
1023 *   const double offset = 2;
1024 *   for (uint i = 0; i < n_elements; ++i)
1025 *   boundary_data[i] =
1026 *   -amplitude * ((i - left_zero) * (i - right_zero)) + offset;
1027 *   }
1028 *  
1029 *  
1030 *   int
1031 *   main()
1032 *   {
1033 *   std::cout << "Boundary participant: starting... \n";
1034 *  
1035 * @endcode
1036 *
1037 * Configuration
1038 *
1039 * @code
1040 *   const std::string configFileName("precice-config.xml");
1041 *   const std::string solverName("boundary-participant");
1042 *   const std::string meshName("boundary-mesh");
1043 *   const std::string dataWriteName("boundary-data");
1044 *  
1045 * @endcode
1046 *
1047 * Adjust to MPI rank and size for parallel computation
1048 *
1049 * @code
1050 *   const int commRank = 0;
1051 *   const int commSize = 1;
1052 *  
1053 *   precice::SolverInterface precice(solverName,
1054 *   configFileName,
1055 *   commRank,
1056 *   commSize);
1057 *  
1058 *   const int meshID = precice.getMeshID(meshName);
1059 *   const int dimensions = precice.getDimensions();
1060 *   const int numberOfVertices = 6;
1061 *  
1062 *   const int dataID = precice.getDataID(dataWriteName, meshID);
1063 *  
1064 * @endcode
1065 *
1066 * Set up data structures
1067 *
1068 * @code
1069 *   std::vector<double> writeData(numberOfVertices);
1070 *   std::vector<int> vertexIDs(numberOfVertices);
1071 *   std::vector<double> vertices(numberOfVertices * dimensions);
1072 *  
1073 * @endcode
1074 *
1075 * Define a boundary mesh
1076 *
1077 * @code
1078 *   std::cout << "Boundary participant: defining boundary mesh \n";
1079 *   const double length = 2;
1080 *   const double xCoord = 1;
1081 *   const double deltaY = length / (numberOfVertices - 1);
1082 *   for (int i = 0; i < numberOfVertices; ++i)
1083 *   for (int j = 0; j < dimensions; ++j)
1084 *   {
1085 *   const unsigned int index = dimensions * i + j;
1086 * @endcode
1087 *
1088 * The x-coordinate is always 1, i.e., the boundary is parallel to the
1089 * y-axis. The y-coordinate is descending from 1 to -1.
1090 *
1091 * @code
1092 *   if (j == 0)
1093 *   vertices[index] = xCoord;
1094 *   else
1095 *   vertices[index] = 1 - deltaY * i;
1096 *   }
1097 *  
1098 * @endcode
1099 *
1100 * Pass the vertices to preCICE
1101 *
1102 * @code
1103 *   precice.setMeshVertices(meshID,
1104 *   numberOfVertices,
1105 *   vertices.data(),
1106 *   vertexIDs.data());
1107 *  
1108 * @endcode
1109 *
1110 * initialize the Solverinterface
1111 *
1112 * @code
1113 *   double dt = precice.initialize();
1114 *  
1115 * @endcode
1116 *
1117 * Start time loop
1118 *
1119 * @code
1120 *   const double end_time = 1;
1121 *   double time = 0;
1122 *   while (precice.isCouplingOngoing())
1123 *   {
1124 * @endcode
1125 *
1126 * Generate new boundary data
1127 *
1128 * @code
1129 *   define_boundary_values(writeData, time, end_time);
1130 *  
1131 *   {
1132 *   std::cout << "Boundary participant: writing coupling data \n";
1133 *   precice.writeBlockScalarData(dataID,
1134 *   numberOfVertices,
1135 *   vertexIDs.data(),
1136 *   writeData.data());
1137 *   }
1138 *  
1139 *   dt = precice.advance(dt);
1140 *   std::cout << "Boundary participant: advancing in time\n";
1141 *  
1142 *   time += dt;
1143 *   }
1144 *  
1145 *   std::cout << "Boundary participant: closing...\n";
1146 *  
1147 *   return 0;
1148 *   }
1149 * @endcode
1150
1151
1152*/
Definition point.h:112
Point< 3 > center
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), 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.
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
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.)
Definition advection.h:75
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:189
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation