From 97badccdc487548e6214746a9be607b6aa319be7 Mon Sep 17 00:00:00 2001 From: kanschat Date: Wed, 8 Aug 2012 16:33:29 +0000 Subject: [PATCH] ip residuals git-svn-id: https://svn.dealii.org/trunk@25780 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/integrators/laplace.h | 201 ++++++++++++++---- 1 file changed, 165 insertions(+), 36 deletions(-) diff --git a/deal.II/include/deal.II/integrators/laplace.h b/deal.II/include/deal.II/integrators/laplace.h index 08e5f1c929..802943abfa 100644 --- a/deal.II/include/deal.II/integrators/laplace.h +++ b/deal.II/include/deal.II/integrators/laplace.h @@ -175,9 +175,8 @@ namespace LocalIntegrators } /** - * Weak boundary condition for the Laplace operator by Nitsche, vector - * valued version, namely on the face F - * the vector + * Weak boundary condition for the Laplace operator by Nitsche, scalar + * version, namely on the face F the vector * @f[ * \int_F \Bigl(\gamma (u-g) v - \partial_n u v - (u-g) \partial_n v\Bigr)\;ds. * @f] @@ -196,40 +195,38 @@ namespace LocalIntegrators void nitsche_residual ( Vector& result, const FEValuesBase& fe, - const VectorSlice > >& input, - const VectorSlice > > >& Dinput, - const VectorSlice > >& data, + const std::vector& input, + const std::vector >& Dinput, + const std::vector& data, double penalty, double factor = 1.) { 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); - AssertVectorVectorDimension(Dinput, n_comp, fe.n_quadrature_points); - AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points); + AssertDimension(input.size(), fe.n_quadrature_points); + AssertDimension(Dinput.size(), fe.n_quadrature_points); + AssertDimension(data.size(), fe.n_quadrature_points); for (unsigned k=0;k& n = fe.normal_vector(k); for (unsigned i=0;iF the vector + * Weak boundary condition for the Laplace operator by Nitsche, vector + * valued version, namely on the face F + * the vector * @f[ * \int_F \Bigl(\gamma (u-g) v - \partial_n u v - (u-g) \partial_n v\Bigr)\;ds. * @f] @@ -248,31 +245,33 @@ namespace LocalIntegrators void nitsche_residual ( Vector& result, const FEValuesBase& fe, - const std::vector& input, - const std::vector >& Dinput, - const std::vector& data, + const VectorSlice > >& input, + const VectorSlice > > >& Dinput, + const VectorSlice > >& data, double penalty, double factor = 1.) { const unsigned int n_dofs = fe.dofs_per_cell; - AssertDimension(input.size(), fe.n_quadrature_points); - AssertDimension(Dinput.size(), fe.n_quadrature_points); - AssertDimension(data.size(), fe.n_quadrature_points); + const unsigned int n_comp = fe.get_fe().n_components(); + AssertVectorVectorDimension(input, n_comp, fe.n_quadrature_points); + AssertVectorVectorDimension(Dinput, n_comp, fe.n_quadrature_points); + AssertVectorVectorDimension(data, n_comp, fe.n_quadrature_points); for (unsigned k=0;k& n = fe.normal_vector(k); for (unsigned i=0;i + void + ip_residual( + Vector& result1, + Vector& result2, + const FEValuesBase& fe1, + const FEValuesBase& fe2, + const std::vector& input1, + const std::vector >& Dinput1, + const std::vector& input2, + const std::vector >& Dinput2, + double pen, + double int_factor = 1., + double ext_factor = -1.) +{ + Assert(fe1.get_fe().n_components() == 1, + ExcDimensionMismatch(fe1.get_fe().n_components(), 1)); + Assert(fe2.get_fe().n_components() == 1, + ExcDimensionMismatch(fe2.get_fe().n_components(), 1)); + + const double nui = int_factor; + const double nue = (ext_factor < 0) ? int_factor : ext_factor; + const double penalty = .5 * pen * (nui + nue); + + const unsigned int n_dofs = fe1.dofs_per_cell; + + for (unsigned k=0;k& n = fe1.normal_vector(k); + + for (unsigned i=0;i& Dvi = fe1.shape_grad(i,k); + const double dnvi = Dvi * n; + const double ve = fe2.shape_value(i,k); + const Tensor<1,dim>& Dve = fe2.shape_grad(i,k); + const double dnve = Dve * n; + + const double ui = input1[k]; + const Tensor<1,dim>& Dui = Dinput1[k]; + const double dnui = Dui * n; + const double ue = input2[k]; + const Tensor<1,dim>& Due = Dinput2[k]; + const double dnue = Due * n; + + result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi); + result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue); + result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve); + result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve); + } + } +} + + +/** + * Vector-valued residual term for the symmetric interior penalty method. + * + * @ingroup Integrators + * @author Guido Kanschat + * @date 2012 + */ + template + 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 n_comp = fe1.get_fe().n_components(); + const unsigned int n1 = fe1.dofs_per_cell; + + AssertVectorVectorDimension(input1, n_comp, fe1.n_quadrature_points); + AssertVectorVectorDimension(Dinput1, n_comp, fe1.n_quadrature_points); + AssertVectorVectorDimension(input2, n_comp, fe2.n_quadrature_points); + AssertVectorVectorDimension(Dinput2, n_comp, fe2.n_quadrature_points); + + const double nui = int_factor; + const double nue = (ext_factor < 0) ? int_factor : ext_factor; + const double penalty = .5 * pen * (nui + nue); + + + for (unsigned k=0;k& n = fe1.normal_vector(k); + + for (unsigned i=0;i& Dvi = fe1.shape_grad_component(i,k,d); + const double dnvi = Dvi * n; + const double ve = fe2.shape_value_component(i,k,d); + const Tensor<1,dim>& Dve = fe2.shape_grad_component(i,k,d); + const double dnve = Dve * n; + + const double ui = input1[d][k]; + const Tensor<1,dim>& Dui = Dinput1[d][k]; + const double dnui = Dui * n; + const double ue = input2[d][k]; + const Tensor<1,dim>& Due = Dinput2[d][k]; + const double dnue = Due * n; + + result1(i) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+penalty*ui*vi); + result1(i) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-penalty*vi*ue); + result2(i) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-penalty*ui*ve); + result2(i) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+penalty*ue*ve); + } + } + } + + /** * Auxiliary function computing the penalty parameter for interior -- 2.39.5