From 6b0cc0df8bbae6950027db61ac5f2c3db11e2169 Mon Sep 17 00:00:00 2001 From: kanschat Date: Mon, 4 Mar 2013 13:13:57 +0000 Subject: [PATCH] add residuals git-svn-id: https://svn.dealii.org/trunk@28721 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/deal.II/integrators/elasticity.h | 181 +++++++++++++++++- 1 file changed, 175 insertions(+), 6 deletions(-) diff --git a/deal.II/include/deal.II/integrators/elasticity.h b/deal.II/include/deal.II/integrators/elasticity.h index 29dfc868c7..8fe0766f8b 100644 --- a/deal.II/include/deal.II/integrators/elasticity.h +++ b/deal.II/include/deal.II/integrators/elasticity.h @@ -1,7 +1,7 @@ //--------------------------------------------------------------------------- // $Id$ // -// Copyright (C) 2010, 2011, 2012 by the deal.II authors +// Copyright (C) 2010, 2011, 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 @@ -35,10 +35,11 @@ namespace LocalIntegrators namespace Elasticity { /** - * Scalar product of symmetric gradients. + * The linear elasticity operator in weak form, namely double + * contraction of symmetric gradients. * * \f[ - * (\varepsilon(u), \varepsilon(v)) + * \int_Z \varepsilon(u): \varepsilon(v)\,dx * \f] */ template @@ -68,9 +69,48 @@ namespace LocalIntegrators /** - * The weak boundary condition - * of Nitsche type for - * symmetric gradients. + * Vector-valued residual operator for linear elasticity in weak form + * + * \f[ + * - \int_Z \varepsilon(u): \varepsilon(v) \,dx + * \f] + */ + template + inline void + cell_residual ( + Vector &result, + const FEValuesBase &fe, + const VectorSlice > > > &input, + double factor = 1.) + { + const unsigned int nq = fe.n_quadrature_points; + const unsigned int n_dofs = fe.dofs_per_cell; + AssertDimension(fe.get_fe().n_components(), dim); + + AssertVectorVectorDimension(input, dim, fe.n_quadrature_points); + Assert(result.size() == n_dofs, ExcDimensionMismatch(result.size(), n_dofs)); + + for (unsigned int k=0; k inline void nitsche_matrix ( @@ -111,6 +151,64 @@ namespace LocalIntegrators } } + /** + * Weak boundary condition for the elasticity operator by Nitsche, + * namely on the face F the vector + * @f[ + * \int_F \Bigl(\gamma (u-g) \cdot v - n^T \epsilon(u) v - (u-g) \epsilon(v) n^T\Bigr)\;ds. + * @f] + * + * Here, u is the finite element function whose values and + * gradient are given in the arguments input and + * Dinput, respectively. g is the inhomogeneous + * boundary value in the argument data. $n$ is the outer + * normal vector and $\gamma$ is the usual penalty parameter. + * + * @author Guido Kanschat + * @date 2013 + */ + template + void nitsche_residual ( + Vector &result, + const FEValuesBase &fe, + const VectorSlice > > &input, + const VectorSlice > > > &Dinput, + const VectorSlice > > &data, + double penalty, + double factor = 1.) + { + const unsigned int n_dofs = fe.dofs_per_cell; + AssertVectorVectorDimension(input, dim, fe.n_quadrature_points); + AssertVectorVectorDimension(Dinput, dim, fe.n_quadrature_points); + AssertVectorVectorDimension(data, dim, fe.n_quadrature_points); + + for (unsigned int k=0; k &n = fe.normal_vector(k); + for (unsigned int i=0; i + void + ip_residual( + Vector &result1, + Vector &result2, + const FEValuesBase &fe1, + const FEValuesBase &fe2, + const VectorSlice > > &input1, + const VectorSlice > > > &Dinput1, + const VectorSlice > > &input2, + const VectorSlice > > > &Dinput2, + double pen, + double int_factor = 1., + double ext_factor = -1.) + { + const unsigned int n1 = fe1.dofs_per_cell; + + AssertDimension(fe1.get_fe().n_components(), dim); + AssertDimension(fe2.get_fe().n_components(), dim); + AssertVectorVectorDimension(input1, dim, fe1.n_quadrature_points); + AssertVectorVectorDimension(Dinput1, dim, fe1.n_quadrature_points); + AssertVectorVectorDimension(input2, dim, fe2.n_quadrature_points); + AssertVectorVectorDimension(Dinput2, dim, fe2.n_quadrature_points); + + const double nu1 = int_factor; + const double nu2 = (ext_factor < 0) ? int_factor : ext_factor; + const double penalty = .5 * pen * (nu1 + nu2); + + + for (unsigned int k=0; k &n = fe1.normal_vector(k); + + for (unsigned int i=0; i