From 946ccaf1cbcd4e5f9ee1d700fa6ae63fc621f705 Mon Sep 17 00:00:00 2001 From: bangerth Date: Tue, 9 Apr 2013 19:48:19 +0000 Subject: [PATCH] Add a 1d mesh worker test. git-svn-id: https://svn.dealii.org/trunk@29230 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/integrators/mesh_worker_1d_dg.cc | 434 ++++++++++++++++++ .../integrators/mesh_worker_1d_dg/cmp/generic | 50 ++ 2 files changed, 484 insertions(+) create mode 100644 tests/integrators/mesh_worker_1d_dg.cc create mode 100644 tests/integrators/mesh_worker_1d_dg/cmp/generic diff --git a/tests/integrators/mesh_worker_1d_dg.cc b/tests/integrators/mesh_worker_1d_dg.cc new file mode 100644 index 0000000000..97feaf0a31 --- /dev/null +++ b/tests/integrators/mesh_worker_1d_dg.cc @@ -0,0 +1,434 @@ +//---------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2000, 2001, 2003, 2004, 2007, 2008, 2009, 2010, 2012, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------------------------------------------------- + +// Test that we can use MeshWorker also in 1d. test by Scott Miller + +#include "../tests.h" +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +using namespace dealii; + + +//! Solve the advection equation: \dot{u} + \div{\mathbf{c} u} = 0.0, +//! with c={1}, {1,0}, {1,0,0} in d=1,2,3 + +//! Domain x \in [0,1], y,z \in [0,0.01] +//! Initial condition: u = 0.0 +//! Boundary condition at x=0: u=1 + +//! Use a DG formulation with upwind fluxes + + +namespace Advection +{ +using namespace dealii; + +/******************************************** + * ADVECTION PROBLEM + ********************************************/ +template +class AdvectionProblem +{ +public: + AdvectionProblem (); + ~AdvectionProblem (); + void run (); + +private: + + const MappingQ1 mapping; + + void setup_system (); + + void integrate_cell_term (MeshWorker::DoFInfo& dinfo, + MeshWorker::IntegrationInfo& info); + + void integrate_boundary_term (MeshWorker::DoFInfo& dinfo, + MeshWorker::IntegrationInfo& info); + + void integrate_face_term (MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + MeshWorker::IntegrationInfo& info1, + MeshWorker::IntegrationInfo& info2); + + void output_results (int timestep) const; + + void create_grid (); + + void assemble_rhs (Vector& solution, + Vector& residual); + + // For problems with non-diagonal mass matrices +// void assemble_mass_matrix_and_multiply (Vector& solution, +// Vector& residual); + + const double wavespeed; + + + // DATA: + Triangulation triangulation; + DoFHandler dof_handler; + + FESystem fe; + + Vector solution, stage; + + const FEValuesExtractors::Scalar upos; +}; + + +template +AdvectionProblem::AdvectionProblem () +: + mapping(), + wavespeed(1.0), + dof_handler (triangulation), + fe (FE_DGQ(0), 1),// p=0, and solving for a scalar + upos(0) +{} + +template +AdvectionProblem::~AdvectionProblem () +{ + dof_handler.clear (); +} + + +template < > +void AdvectionProblem<1>::create_grid() +{ + double ll_x=0.; + double ur_x=1.; + + int n_cells_x = 10; + + const Point<1> LowerLeft (ll_x), + UpperRight (ur_x); + + // Define the subdivisions in the x1 and x2 coordinates. + std::vector subdivisions(1); + subdivisions[0] = n_cells_x; + + GridGenerator::subdivided_hyper_rectangle(triangulation, + subdivisions, + LowerLeft, + UpperRight, + true); + +}//create_grid() + +template < > +void AdvectionProblem<2>::create_grid() +{ + // dim==1 implemented: + double ll_x=0., ll_y=0.; + double ur_x=1., ur_y=0.01; + + int n_cells_x = 100; + int n_cells_y = 1; + + const Point<2> LowerLeft (ll_x, ll_y), + UpperRight (ur_x, ur_y ); + + // Define the subdivisions in the x1 and x2 coordinates. + std::vector subdivisions(2); + subdivisions[0] = n_cells_x; + subdivisions[1] = n_cells_y; + + GridGenerator::subdivided_hyper_rectangle(triangulation, + subdivisions, + LowerLeft, + UpperRight, + true); + +}//create_grid() + +template <> +void AdvectionProblem<3>::create_grid() +{ + double ll_x=0., ll_y=0., ll_z=0.; + double ur_x=1., ur_y=0.01, ur_z=0.01; + + int n_cells_x = 100; + int n_cells_y = 1; + int n_cells_z = 1; + + const Point<3> LowerLeft (ll_x, ll_y, ll_z), + UpperRight (ur_x, ur_y, ur_z); + + // Define the subdivisions in the x1 and x2 coordinates. + std::vector subdivisions(3); + subdivisions[0] = n_cells_x; + subdivisions[1] = n_cells_y; + subdivisions[2] = n_cells_z; + + GridGenerator::subdivided_hyper_rectangle(triangulation, + subdivisions, + LowerLeft, + UpperRight, + true); + +}//create_grid() + +template +void AdvectionProblem::setup_system () +{ + dof_handler.distribute_dofs (fe); + solution.reinit(dof_handler.n_dofs()); + stage.reinit(dof_handler.n_dofs()); + + solution = 0.0; + stage = 0.0; +}//setup_system + + +template +void AdvectionProblem::assemble_rhs (Vector &solution, + Vector& residual) +{ + const unsigned int n_gauss_points = std::ceil(((2.0*fe.degree) +1)/2); + + MeshWorker::IntegrationInfoBox info_box; + + info_box.initialize_gauss_quadrature(n_gauss_points, + n_gauss_points, + n_gauss_points); + + info_box.initialize_update_flags(); + UpdateFlags update_flags = update_quadrature_points | + update_values | + update_gradients; + + info_box.add_update_flags(update_flags, true, true, true, true); + + NamedData* > solution_data; + + Vector* u = &solution; + + solution_data.add(u, "solution"); + info_box.cell_selector.add("solution", true, true, false); + info_box.boundary_selector.add("solution", true, false, false); + info_box.face_selector.add("solution", true, false, false); + + info_box.initialize(fe, mapping, solution_data); + +//deallog<<"\nWe are now going to attend construction of MeshWorker::DoFInfo..."< dof_info(dof_handler); +//deallog<<"\nApparently it DoFInfo was constructed fine!"< > assembler; + NamedData* > data; + Vector* rhs = &residual; + data.add(rhs, "Residual"); + assembler.initialize(data); + + MeshWorker::loop, MeshWorker::IntegrationInfoBox > + (dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + std_cxx1x::bind(&AdvectionProblem::integrate_cell_term, + this, std_cxx1x::_1, std_cxx1x::_2), + std_cxx1x::bind(&AdvectionProblem::integrate_boundary_term, + this, std_cxx1x::_1, std_cxx1x::_2), + std_cxx1x::bind(&AdvectionProblem::integrate_face_term, + this, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3, std_cxx1x::_4), + assembler, true); + +}//assemble_system + +template +void AdvectionProblem::integrate_cell_term (MeshWorker::DoFInfo& dinfo, + MeshWorker::IntegrationInfo& info) +{ + const FEValuesBase& fe_v = info.fe_values(); + + Vector& cell_rhs = dinfo.vector(0).block(0); + + const unsigned int dofs_per_cell = fe_v.dofs_per_cell; + const unsigned int n_q_points = fe_v.n_quadrature_points; + + FullMatrix cell_matrix(dofs_per_cell,dofs_per_cell); + cell_matrix = 0.0; + + const std::vector > &values = info.values[0]; + + std::vector u(values[0]); + + for (unsigned int q_point=0; q_point +void AdvectionProblem::integrate_boundary_term (MeshWorker::DoFInfo& dinfo, + MeshWorker::IntegrationInfo& info) +{ + const unsigned int boundary_id = dinfo.face->boundary_indicator(); + + // We only have a non-zero boundary contribution at the + // x=0 boundary + if(boundary_id != 0) + return; + + const FEValuesBase& fe_v = info.fe_values(); + + Vector& cell_rhs = dinfo.vector(0).block(0); + + const unsigned int dofs_per_cell = fe_v.dofs_per_cell; + const unsigned int n_q_points = fe_v.n_quadrature_points; + + double boundary_flux = -1.0; + + for (unsigned int q_point=0; q_point +void AdvectionProblem::integrate_face_term (MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + MeshWorker::IntegrationInfo& info1, + MeshWorker::IntegrationInfo& info2) +{ + const FEValuesBase& fe_v_1 = info1.fe_values(); + const FEValuesBase& fe_v_2 = info2.fe_values(); + + Vector& cell_vector_1 = dinfo1.vector(0).block(0); + Vector& cell_vector_2 = dinfo2.vector(0).block(0); + + const unsigned int dofs_per_cell = fe_v_1.dofs_per_cell; + const unsigned int n_q_points = fe_v_1.n_quadrature_points; + + std::vector &u_1 = info1.values[0][0]; + std::vector &u_2 = info1.values[0][0]; + + double flux; + + for (unsigned int q_point=0; q_point0) + flux = u_1[q_point]*fe_v_1.normal_vector(q_point)[0]; + else + flux = u_2[q_point]*fe_v_1.normal_vector(q_point)[0]; + + for (unsigned int i=0; i +void AdvectionProblem::output_results (int timestep) const +{ + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + + std::vector solution_names(1, "u"); + + std::vector + interpretation (1, DataComponentInterpretation::component_is_scalar); + + data_out.add_data_vector (solution, + solution_names, + DataOut >::type_automatic, + interpretation); + + data_out.build_patches (fe.degree); + data_out.write_gnuplot (deallog.get_file_stream()); + +}//output_results + + +template +void AdvectionProblem::run () +{ + // Make the mesh + create_grid(); + + // Setup the system + setup_system(); + + deallog << "\tNumber of active cells: " + << triangulation.n_active_cells() << std::endl; + + deallog << "\tNumber of degrees of freedom: " + << dof_handler.n_dofs() << std::endl; + + const double delta_t = 0.005; + const double n_dt = 10; + + double inv_cell_vol = 1.0/std::pow(0.01, dim); + + for (unsigned int dt=0; dt advection_problem; + advection_problem.run (); +} diff --git a/tests/integrators/mesh_worker_1d_dg/cmp/generic b/tests/integrators/mesh_worker_1d_dg/cmp/generic new file mode 100644 index 0000000000..65c3e36044 --- /dev/null +++ b/tests/integrators/mesh_worker_1d_dg/cmp/generic @@ -0,0 +1,50 @@ + +DEAL:: Number of active cells: 10 +DEAL:: Number of degrees of freedom: 10 +# This file was generated by the deal.II library. + + +# +# For a description of the GNUPLOT format see the GNUPLOT manual. +# +# +0.00000 0.999023 +0.100000 0.999023 + + +0.100000 0.989258 +0.200000 0.989258 + + +0.200000 0.945312 +0.300000 0.945312 + + +0.300000 0.828125 +0.400000 0.828125 + + +0.400000 0.623047 +0.500000 0.623047 + + +0.500000 0.376953 +0.600000 0.376953 + + +0.600000 0.171875 +0.700000 0.171875 + + +0.700000 0.0546875 +0.800000 0.0546875 + + +0.800000 0.0107422 +0.900000 0.0107422 + + +0.900000 0.000976562 +1.00000 0.000976562 + + -- 2.39.5