From: hartmann Date: Fri, 9 Nov 2001 13:14:19 +0000 (+0000) Subject: Linear transport in non-conservative form. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fbb6e29445f52d5fc5fcc7fd87914790aa9fa456;p=dealii-svn.git Linear transport in non-conservative form. git-svn-id: https://svn.dealii.org/trunk@5187 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc new file mode 100644 index 0000000000..ff3e5eeaf1 --- /dev/null +++ b/deal.II/examples/step-12/step-12.cc @@ -0,0 +1,716 @@ +/* $Id$ */ +/* Author: Ralf Hartmann, University of Heidelberg, 2000 */ + + // The first few files have already + // been covered in previous examples + // and will thus not be further + // commented on. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + + +#include + + +template +class Beta +{ + public: + Beta () {}; + + void value_list (const std::vector > &points, + std::vector > &values) const; +}; + + +template +class RHS: public Function +{ + public: + RHS() {}; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component=0) const; +}; + + +template +class BoundaryFunction: public Function +{ + public: + BoundaryFunction() {}; + + virtual void value_list (const std::vector > &points, + std::vector &values, + const unsigned int component=0) const; +}; + +template +class DGAssembler +{ + public: + DGAssembler() {}; + + void assemble_cell_term(const FEValuesBase& fe_v, + FullMatrix &cell_matrix, + Vector &cell_vector); + + void assemble_face_term(const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &cell_matrix, + FullMatrix &cell_inflow_matrix, + Vector &cell_vector); + + private: + Beta beta_function; + RHS rhs_function; + BoundaryFunction boundary_function; +}; + + // The main class is again almost + // unchanged. Two additions, however, + // are made: we have added the + // ``refine'' function, which is used + // to adaptively refine the grid + // (instead of the global refinement + // in the previous examples), and a + // variable which will hold the + // constraints associated to the + // hanging nodes. +template +class TransportProblem +{ + public: + TransportProblem (); + ~TransportProblem (); + + void run (); + + private: + void setup_system (); + void assemble_system (); + void solve (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + MappingQ1 mapping; + + // We need a finite element + // again. This time, we will want + // to use quadratic polynomials + // (but this is only specified in + // the constructor): + FE_DGQ fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector right_hand_side; + + DGAssembler dg_assembler; +}; + + +template <> +void Beta<2>::value_list(const std::vector > &points, + std::vector > &values) const +{ + Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size())); + for (unsigned int i=0; i &p=points[i]; + Point<2> &beta=values[i]; + + beta(0)=-p(1); + beta(1)=p(0); + beta/=sqrt(beta*beta); + } +} + + + + +template +void RHS::value_list(const std::vector > &, + std::vector &values, + const unsigned int) const +{ + for (unsigned int i=0; i +void BoundaryFunction::value_list(const std::vector > &points, + std::vector &values, + const unsigned int) const +{ + Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size())); + for (unsigned int i=0; i +void DGAssembler::assemble_cell_term(const FEValuesBase& fe_v, + FullMatrix &cell_matrix, + Vector &cell_vector) +{ + const vector > > &grad_v = fe_v.get_shape_grads (); + const FullMatrix &v = fe_v.get_shape_values (); + const vector &JxW = fe_v.get_JxW_values (); + + vector > beta (fe_v.n_quadrature_points); + vector rhs (fe_v.n_quadrature_points); + + beta_function.value_list (fe_v.get_quadrature_points(), beta); + rhs_function.value_list (fe_v.get_quadrature_points(), rhs); + + for (unsigned int point=0; point +void DGAssembler::assemble_face_term(const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &cell_matrix, + FullMatrix &cell_inflow_matrix, + Vector &cell_vector) +{ + DoFHandler::face_iterator face=fe_v.get_face(); + + const FullMatrix &v = fe_v.get_shape_values (); + const FullMatrix &v_neighbor = fe_v_neighbor.get_shape_values (); + const vector &JxW = fe_v.get_JxW_values (); + const vector > &normals = fe_v.get_normal_vectors (); + + vector > beta (fe_v.n_quadrature_points); + vector g(fe_v.n_quadrature_points); + + beta_function.value_list (fe_v.get_quadrature_points(), beta); + + if (face->at_boundary()) + boundary_function.value_list (fe_v.get_quadrature_points(), g); + + for (unsigned int point=0; pointat_boundary()) + for (unsigned int k=0; k +TransportProblem::TransportProblem () : + fe (1), + dof_handler (triangulation) +{} + + +template +TransportProblem::~TransportProblem () +{ + dof_handler.clear (); +}; + + + +template +void TransportProblem::setup_system () +{ + // To distribute degrees of + // freedom, the ``dof_handler'' + // variable takes only the finite + // element object. In this case, it + // will distribute four degrees of + // freedom per cell. + dof_handler.distribute_dofs (fe); + + sparsity_pattern.reinit (dof_handler.n_dofs(), + dof_handler.n_dofs(), + dof_handler.max_couplings_between_dofs()); + DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern); + sparsity_pattern.compress(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + right_hand_side.reinit (dof_handler.n_dofs()); +}; + + + +template +void TransportProblem::assemble_system () +{ + // See Cockburn paper for the proper quadrature. + QGauss2 quadrature; + QGauss2 face_quadrature; + + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + vector dofs (dofs_per_cell); + vector dofs_neighbor (dofs_per_cell); + + UpdateFlags update_flags = UpdateFlags(update_values + | update_gradients + | update_q_points + | update_JxW_values); + + UpdateFlags face_update_flags = UpdateFlags(update_values + | update_q_points + | update_JxW_values + | update_normal_vectors); + + + FEValues fe_v ( + mapping, fe, quadrature, update_flags); + FEFaceValues fe_v_face ( + mapping, fe, face_quadrature, face_update_flags); + FESubfaceValues fe_v_subface ( + mapping, fe, face_quadrature, face_update_flags); + FEFaceValues fe_v_face_neighbor ( + mapping, fe, face_quadrature, UpdateFlags(update_values | update_default)); + FESubfaceValues fe_v_subface_neighbor ( + mapping, fe, face_quadrature, UpdateFlags(update_values | update_default)); + + // includes the u and v terms + FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); + // includes u_hat and v terms + FullMatrix cell_inflow_matrix (dofs_per_cell, dofs_per_cell); + + Vector cell_vector (dofs_per_cell); + + DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + DoFHandler::face_iterator face; + DoFHandler::cell_iterator neighbor; + DoFHandler::cell_iterator neighbor_child; + + for (;cell!=endc; ++cell) + { + // re-init fe values for this cell + fe_v.reinit (cell); + + cell_matrix.clear (); + cell_vector.clear (); + + dg_assembler.assemble_cell_term(fe_v, + cell_matrix, + cell_vector); + + cell->get_dof_indices (dofs); + + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + face = cell->face(face_no); + + cell_inflow_matrix.clear(); + + if (face->at_boundary()) + { + fe_v_face.reinit (cell, face_no); + + dg_assembler.assemble_face_term(fe_v_face, + fe_v_face, + cell_matrix, + cell_inflow_matrix, + cell_vector); + } + else // if (!face->at_boundary()) + { + Assert (cell->neighbor(face_no).state() == valid, ExcInternalError()); + neighbor = cell->neighbor(face_no); + + if (face->has_children()) // i.e. neighbor is one level more refined than cell + { + // store which number #cell# has in the + // list of neighbors of #neighbor# + const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); + + + // loop over all subfaces + for (unsigned int subface_no=0; subface_no::subfaces_per_face; + ++subface_no) + { + // get an iterator pointing to the + // cell behind the present subface + neighbor_child = neighbor->child(GeometryInfo:: + child_cell_on_face(neighbor2,subface_no)); + Assert (neighbor_child->face(neighbor2) == face->child(subface_no), + ExcInternalError()); + Assert (!neighbor_child->has_children(), ExcInternalError()); + + fe_v_subface.reinit (cell, face_no, subface_no); + fe_v_face_neighbor.reinit (neighbor_child, neighbor2); + + cell_inflow_matrix.clear(); + + dg_assembler.assemble_face_term(fe_v_subface, + fe_v_face_neighbor, + cell_matrix, + cell_inflow_matrix, + cell_vector); + + // get indices of dofs of neighbor_child cell + neighbor_child->get_dof_indices (dofs_neighbor); + + // distribute cell matrix + for (unsigned int i=0; ihas_children()) + { + if (neighbor->level() == cell->level()) + { + // store which number #cell# has in the + // list of neighbors of #neighbor# + const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); + + fe_v_face.reinit (cell, face_no); + fe_v_face_neighbor.reinit (neighbor, neighbor2); + + dg_assembler.assemble_face_term(fe_v_face, + fe_v_face_neighbor, + cell_matrix, + cell_inflow_matrix, + cell_vector); + } + else // if (neighbor->level() < cell->level()) i.e. neighbor is one level coarser than cell + { + Assert(neighbor->level() < cell->level(), ExcInternalError()); + + const std::pair faceno_subfaceno= + cell->neighbor_of_coarser_neighbor(face_no); + const unsigned int neighbor_face_no=faceno_subfaceno.first, + neighbor_subface_no=faceno_subfaceno.second; + + Assert (neighbor->neighbor(neighbor_face_no) + ->child(GeometryInfo::child_cell_on_face( + face_no,neighbor_subface_no)) == cell, ExcInternalError()); + + // now 'neighbor_face_no' stores the number + // of a face in the list of faces of 'neighbor'. + // This face has got a subface that is + // between 'cell' and 'neighbor'. + // 'neighbor_subface_no' stores the number + // of this subface in the list of subfaces of this + // face 'neighbor->face(neighbor_face_no)' + // that is between 'cell' and 'neighbor' + fe_v_face.reinit (cell, face_no); + fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no, + neighbor_subface_no); + + dg_assembler.assemble_face_term(fe_v_face, + fe_v_subface_neighbor, + cell_matrix, + cell_inflow_matrix, + cell_vector); + } // else // if (neighbor->level() < cell->level()) + + // get indices of dofs of neighbor_child cell + neighbor->get_dof_indices (dofs_neighbor); + + // distribute cell_inflow_matrix + for (unsigned int i=0; ihas_children()) + } // else // if (!face->at_boundary()) + } //for (face_no...) + + // distribute cell matrix + for (unsigned int i=0; i +void TransportProblem::solve () +{ + SolverControl solver_control (1000, 1e-12); + PrimitiveVectorMemory<> vector_memory; + SolverRichardson<> solver (solver_control, vector_memory); + + PreconditionBlockSSOR preconditioner; + preconditioner.initialize(system_matrix, fe.dofs_per_cell); + preconditioner.invert_diagblocks(); + + solver.solve (system_matrix, solution, right_hand_side, + preconditioner); +}; + + +template +void TransportProblem::refine_grid () +{ + Vector estimated_error_per_cell (triangulation.n_active_cells()); + + FunctionMap::type neumann_boundary; + + KellyErrorEstimator::estimate (dof_handler, + QGauss3(), + neumann_boundary, + solution, + estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + + triangulation.execute_coarsening_and_refinement (); +}; + + + +template +void TransportProblem::output_results (const unsigned int cycle) const +{ + // We want to write the grid in + // each cycle. Here is another way + // to quickly produce a filename + // based on the cycle number. It + // assumes that the numbers `0' + // through `9' are represented + // consecutively in the character + // set (which is the case in all + // known character sets). However, + // this will only work if the cycle + // number is less than ten, which + // we check by an assertion. + std::string filename = "grid-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".eps"; + std::ofstream eps_output (filename.c_str()); + + // Using this filename, we write + // each grid as a postscript file. + GridOut grid_out; + grid_out.write_eps (triangulation, eps_output); + + // output of the solution + filename = "sol-"; + filename += ('0' + cycle); + Assert (cycle < 10, ExcInternalError()); + + filename += ".gnuplot"; + std::ofstream gnuplot_output (filename.c_str()); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "u"); + + data_out.build_patches (); + + data_out.write_gnuplot(gnuplot_output); +}; + + + +template +void TransportProblem::run () +{ + for (unsigned int cycle=0; cycle<3; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + GridGenerator::hyper_cube (triangulation); + + triangulation.refine_global (3); + } + else + // In case this is not the + // first cycle, we want to + // refine the grid. Unlike + // the global refinement + // employed in the last + // example, we now use the + // adaptive procedure + // described in the function + // which we now call: + { + refine_grid (); + }; + + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + + setup_system (); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + assemble_system (); + solve (); + output_results (cycle); + } +} + +int main () +{ + + // The general idea behind the + // layout of this function is as + // follows: let's try to run the + // program as we did before... + try + { + TransportProblem<2> Transport_problem_2d; + Transport_problem_2d.run (); + } + // ...and if this should fail, try + // to gather as much information as + // possible. Specifically, if the + // exception that was thrown is an + // object of a class that is + // derived from the C++ standard + // class ``exception'', then we can + // use the ``what'' member function + // to get a string which describes + // the reason why the exception was + // thrown. + // + // The deal.II exception classes + // are all derived from the + // standard class, and in + // particular, the ``exc.what()'' + // function will return + // approximately the same string as + // would be generated if the + // exception was thrown using the + // ``Assert'' macro. You have seen + // the output of such an exception + // in the previous example, and you + // then know that it contains the + // file and line number of where + // the exception occured, and some + // other information. This is also + // what would be printed in the + // following. + 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; + // We can't do much more than + // printing as much information + // as we can get to, so abort + // with error: + return 1; + } + // If the exception that was thrown + // somewhere was not an object of a + // class derived from the standard + // ``exception'' class, then we + // can't do anything at all. We + // then simply print an error + // message and exit. + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; + + // If we got to this point, there + // was no exception which + // propagated up to the main + // function (maybe there were some, + // but they were caught somewhere + // in the program or the + // library). Therefore, the program + // performed as was expected and we + // can return without error. + return 0; +};