From: kanschat Date: Sun, 3 Jan 2010 15:50:43 +0000 (+0000) Subject: simplify example X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a51136adecb7457f5790a0300b75e4941f8f7886;p=dealii-svn.git simplify example git-svn-id: https://svn.dealii.org/trunk@20274 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-38/doc/intro.dox b/deal.II/examples/step-38/doc/intro.dox index 30db93f368..fdc96eed63 100644 --- a/deal.II/examples/step-38/doc/intro.dox +++ b/deal.II/examples/step-38/doc/intro.dox @@ -5,9 +5,11 @@ This example is devoted to the MeshWorker framework and the discontinuous Galerkin method, or in short: DG method. It -solves the same problem as @ref step_12 "step-12" (see there for a description of the -problem and discretization), but here we use the MeshWorker framework -in order to save reprogramming the cell/face loops. +solves the same problem as @ref step_12 "step-12" (see there for a +description of the problem and discretization), but here we use the +MeshWorker framework in order to save reprogramming the cell/face +loops. We have tried to strip this example of peripheral information +such that the structure becomes more clear. In particular the loops of DG methods turn out to be complex, because for the face terms, we have to distinguish the cases of boundary, @@ -16,3 +18,8 @@ respectively. The MeshWorker framework implements the standard loop over all cells and faces in MeshWorker::loop() and takes care of distinguishing between all the different faces. +There are two things left to do if you use MeshWorker: first, you need +to write the local integrators for your problem. Second, you select +classes from the MeshWorker namespace and combine them to achieve your +goal. + diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc index 122db9727f..090e4edb86 100644 --- a/deal.II/examples/step-38/step-38.cc +++ b/deal.II/examples/step-38/step-38.cc @@ -3,7 +3,7 @@ /* $Id$ */ /* */ -/* Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009 by the deal.II authors */ +/* Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -11,7 +11,7 @@ /* further information on this license. */ // The first few files have already - // been covered in example 12 + // been covered in step-12 // and will thus not be further // commented on. #include @@ -50,27 +50,6 @@ using namespace dealii; // @sect3{Equation data} // - // First we define the classes - // representing the equation-specific - // functions. Both classes, RHS - // and BoundaryValues, are - // derived from the Function - // class. Only the value_list - // function are implemented because - // only lists of function values are - // computed rather than single - // values. -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 BoundaryValues: public Function { @@ -81,80 +60,7 @@ class BoundaryValues: public Function const unsigned int component=0) const; }; - - // The class Beta represents the - // vector valued flow field of the - // linear transport equation and is - // not derived from the Function - // class as we prefer to get function - // values of type Point rather - // than of type - // Vector@. This, because - // there exist scalar products - // between Point and Point as - // well as between Point and - // Tensor, simplifying terms like - // $\beta\cdot n$ and - // $\beta\cdot\nabla v$. - // - // An unnecessary empty constructor - // is added to the class to work - // around a bug in Compaq's cxx - // compiler which otherwise reports - // an error about an omitted - // initializer for an object of - // this class further down. -template -class Beta -{ - public: - Beta () {} - void value_list (const std::vector > &points, - std::vector > &values) const; -}; - - - // The implementation of the - // value_list functions of these - // classes are rather simple. For - // simplicity the right hand side is - // set to be zero but will be - // assembled anyway. -template -void RHS::value_list(const std::vector > &points, - std::vector &values, - const unsigned int) const -{ - // Usually we check whether input - // parameter have the right sizes. - Assert(values.size()==points.size(), - ExcDimensionMismatch(values.size(),points.size())); - - for (unsigned int i=0; i -void Beta::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 beta_function; - const RHS rhs_function; - const BoundaryValues boundary_function; + BoundaryValues boundary_function; }; // @sect4{The local integrators} @@ -258,27 +159,35 @@ void DGIntegrator::cell(CellInfo& info) const const FEValuesBase& fe_v = info.fe(); FullMatrix& local_matrix = info.M1[0].matrix; Vector& local_vector = info.R[0].block(0); - - // With these objects, we continue - // local integration like in step-12. const std::vector &JxW = fe_v.get_JxW_values (); - std::vector > beta (fe_v.n_quadrature_points); - std::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); - + // With these objects, we continue + // local integration like + // always. First, we loop over the + // quadrature points and compute + // the advection vector in the + // current point. for (unsigned int point=0; point beta; + beta(0) = -fe_v.quadrature_point(point)(1); + beta(1) = fe_v.quadrature_point(point)(0); + beta /= beta.norm(); + + for (unsigned int i=0; i::bdry(FaceInfo& info) const const std::vector &JxW = fe_v.get_JxW_values (); const std::vector > &normals = fe_v.get_normal_vectors (); - std::vector > beta (fe_v.n_quadrature_points); std::vector g(fe_v.n_quadrature_points); - beta_function.value_list (fe_v.get_quadrature_points(), beta); boundary_function.value_list (fe_v.get_quadrature_points(), g); for (unsigned int point=0; point beta; + beta(0) = -fe_v.quadrature_point(point)(1); + beta(1) = fe_v.quadrature_point(point)(0); + beta /= beta.norm(); + + const double beta_n=beta * normals[point]; if (beta_n>0) for (unsigned int i=0; i::face(FaceInfo& info1, FaceInfo& info2) const const std::vector &JxW = fe_v.get_JxW_values (); const std::vector > &normals = fe_v.get_normal_vectors (); - - std::vector > beta (fe_v.n_quadrature_points); - beta_function.value_list (fe_v.get_quadrature_points(), beta); - for (unsigned int point=0; point beta; + beta(0) = -fe_v.quadrature_point(point)(1); + beta(1) = fe_v.quadrature_point(point)(0); + beta /= beta.norm(); + + const double beta_n=beta * normals[point]; if (beta_n>0) { // This term we've already