From: kanschat Date: Tue, 4 Sep 2012 08:42:26 +0000 (+0000) Subject: Implement advection operators in weak form X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b7da05410aae61033e6cd80ed7542a5f6c78b8e7;p=dealii-svn.git Implement advection operators in weak form git-svn-id: https://svn.dealii.org/trunk@26230 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/integrators/advection.h b/deal.II/include/deal.II/integrators/advection.h new file mode 100644 index 0000000000..56cdbbcdd4 --- /dev/null +++ b/deal.II/include/deal.II/integrators/advection.h @@ -0,0 +1,390 @@ +//--------------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2010, 2011, 2012 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. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__integrators_advection_h +#define __deal2__integrators_advection_h + + +#include +#include +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace LocalIntegrators +{ +/** + * @brief Local integrators related to advection along a vector field and its DG formulations + * + * All advection operators depend on an advection velocity denoted by + * w in the formulas below. It is denoted as velocity + * in the parameter lists. + * + * The functions cell_matrix() and both upwind_value_matrix() are + * taking the equation in weak form, that is, the directional + * derivative is on the test function. + * + * @ingroup Integrators + * @author Guido Kanschat + * @date 2012 + */ + namespace Advection + { +/** + * Advection along the direction w in weak form + * with derivative on the test function + * \f[ + * m_{ij} = \int_Z u_j\,(\mathbf w \cdot \nabla) v_i \, dx. + * \f] + * + * The FiniteElement in fe may be scalar or vector valued. In + * the latter case, the advection operator is applied to each component + * separately. + * + * @param M: The advection matrix obtained as result + * @param fe: The FEValues object describing the local trial function + * space. #update_values and #update_gradients, and #update_JxW_values + * must be set. + * @param fetest: The FEValues object describing the local test + * function space. #update_values and #update_gradients must be set. + * @param velocity: The advection velocity, a vector of dimension + * dim. Each component may either contain a vector of length + * one, in which case a constant velocity is assumed, or a vector with + * as many entries as quadrature points if the velocity is not constant. + * @param factor is an optional multiplication factor for the result. + * + * @ingroup Integrators + * @author Guido Kanschat + * @date 2012 + */ + template + void cell_matrix ( + FullMatrix& M, + const FEValuesBase& fe, + const FEValuesBase& fetest, + const VectorSlice > >& velocity, + const double factor = 1.) + { + const unsigned int n_dofs = fe.dofs_per_cell; + const unsigned int t_dofs = fetest.dofs_per_cell; + const unsigned int n_components = fe.get_fe().n_components(); + + AssertDimension(velocity.size(), dim); + // If the size of the + // velocity vectors is one, + // then do not increment + // between quadrature points. + const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1; + + if (v_increment == 1) + { + AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points); + } + + AssertDimension(M.n(), n_dofs); + AssertDimension(M.m(), t_dofs); + + for (unsigned k=0;k + inline void + cell_residual ( + Vector& result, + const FEValuesBase& fe, + const std::vector >& input, + const VectorSlice > >& velocity, + double factor = 1.) + { + const unsigned int nq = fe.n_quadrature_points; + const unsigned int n_dofs = fe.dofs_per_cell; + Assert(input.size() == nq, ExcDimensionMismatch(input.size(), nq)); + Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs)); + + AssertDimension(velocity.size(), dim); + const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1; + if (v_increment == 1) + { + AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points); + } + + for (unsigned k=0;k + inline void + cell_residual ( + Vector& result, + const FEValuesBase& fe, + const VectorSlice > > >& input, + const VectorSlice > >& velocity, + double factor = 1.) + { + const unsigned int nq = fe.n_quadrature_points; + const unsigned int n_dofs = fe.dofs_per_cell; + const unsigned int n_comp = fe.get_fe().n_components(); + + AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points); + Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs)); + + AssertDimension(velocity.size(), dim); + const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1; + if (v_increment == 1) + { + AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points); + } + + for (unsigned k=0;kvelocity is + * provided as a VectorSlice, + * having dim vectors, + * one for each velocity + * component. Each of the + * vectors must either have only + * a single entry, if t he + * advection velocity is + * constant, or have an entry + * for each quadrature point. + * + * The finite element can have + * several components, in which + * case each component is + * advected by the same velocity. + */ + template + void upwind_value_matrix( + FullMatrix& M, + const FEValuesBase& fe, + const FEValuesBase& fetest, + const VectorSlice > >& velocity, + double factor = 1.) + { + const unsigned int n_dofs = fe.dofs_per_cell; + const unsigned int t_dofs = fetest.dofs_per_cell; + unsigned int n_components = fe.get_fe().n_components(); + AssertDimension (M.m(), n_dofs); + AssertDimension (M.n(), n_dofs); + + AssertDimension(velocity.size(), dim); + const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1; + if (v_increment == 1) + { + AssertVectorVectorDimension(velocity, dim, fe.n_quadrature_points); + } + + for (unsigned k=0;k 0) + { + for (unsigned i=0;ivelocity is + * provided as a VectorSlice, + * having dim vectors, + * one for each velocity + * component. Each of the + * vectors must either have only + * a single entry, if t he + * advection velocity is + * constant, or have an entry + * for each quadrature point. + * + * The finite element can have + * several components, in which + * case each component is + * advected the same way. + */ + template + void upwind_value_matrix ( + FullMatrix& M11, + FullMatrix& M12, + FullMatrix& M21, + FullMatrix& M22, + const FEValuesBase& fe1, + const FEValuesBase& fe2, + const FEValuesBase& fetest1, + const FEValuesBase& fetest2, + const VectorSlice > >& velocity, + const double factor = 1.) + { + const unsigned int n1 = fe1.dofs_per_cell; + // Multiply the quadrature point + // index below with this factor to + // have simpler data for constant + // velocities. + AssertDimension(velocity.size(), dim); + const unsigned int v_increment = (velocity[0].size() == 1) ? 0 : 1; + if (v_increment == 1) + { + AssertVectorVectorDimension(velocity, dim, fe1.n_quadrature_points); + } + + for (unsigned k=0;k 0) + { + M11(i,j) += dx_nbeta + * fe1.shape_value(j,k) + * fetest1.shape_value(i,k); + M21(i,j) -= dx_nbeta + * fe1.shape_value(j,k) + * fetest2.shape_value(i,k); + } + else + { + M22(i,j) -= dx_nbeta + * fe2.shape_value(j,k) + * fetest2.shape_value(i,k); + M12(i,j) += dx_nbeta + * fe2.shape_value(j,k) + * fetest1.shape_value(i,k); + } + } + else + { + for (unsigned int d=0;d 0) + { + M11(i,j) += dx_nbeta + * fe1.shape_value_component(j,k,d) + * fetest1.shape_value_component(i,k,d); + M21(i,j) -= dx_nbeta + * fe1.shape_value_component(j,k,d) + * fetest2.shape_value_component(i,k,d); + } + else + { + M22(i,j) -= dx_nbeta + * fe2.shape_value_component(j,k,d) + * fetest2.shape_value_component(i,k,d); + M12(i,j) += dx_nbeta + * fe2.shape_value_component(j,k,d) + * fetest1.shape_value_component(i,k,d); + } + } + } + } + } +} + + +DEAL_II_NAMESPACE_CLOSE + +#endif