From 82cd646683574c67b66357b5c2963b2a363fb697 Mon Sep 17 00:00:00 2001 From: kanschat Date: Mon, 18 Oct 2010 23:19:57 +0000 Subject: [PATCH] add Nitsche and IP git-svn-id: https://svn.dealii.org/trunk@22380 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/integrators/laplace.h | 213 ++++++++++++++++-- 1 file changed, 195 insertions(+), 18 deletions(-) diff --git a/deal.II/deal.II/include/integrators/laplace.h b/deal.II/deal.II/include/integrators/laplace.h index 0b05209326..1b26578d5a 100644 --- a/deal.II/deal.II/include/integrators/laplace.h +++ b/deal.II/deal.II/include/integrators/laplace.h @@ -20,6 +20,7 @@ #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -47,10 +48,10 @@ namespace LocalIntegrators * @date 2008, 2009, 2010 */ template - void cell_matrix ( - FullMatrix& M, - const FEValuesBase& fe, - const double factor = 1.) + void cell_matrix ( + FullMatrix& M, + const FEValuesBase& fe, + const double factor = 1.) { const unsigned int n_dofs = fe.dofs_per_cell; const unsigned int n_components = fe.get_fe().n_components(); @@ -63,45 +64,221 @@ namespace LocalIntegrators for (unsigned j=0;jF the matrix + * Weak boundary condition of Nitsche type for the Laplacian, namely on the face F the matrix * @f[ * \int_F \Bigl(\gamma u v - \partial_n u v - u \partial_n v\Bigr)\;ds. * @f] * * Here, γ is the penalty parameter suitably computed * with compute_penalty(). + * + * @author Guido Kanschat + * @date 2008, 2009, 2010 + */ + template + void nitsche_matrix ( + FullMatrix& M, + const FEValuesBase& fe, + double penalty, + double factor = 1.) + { + const unsigned int n_dofs = fe.dofs_per_cell; + const unsigned int n_comp = fe.get_fe().n_components(); + + Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs)); + Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs)); + + for (unsigned k=0;k& n = fe.normal_vector(k); + for (unsigned i=0;iF + * the vector + * @f[ + * \int_F \Bigl(\gamma (u-g) v - \partial_n u v - (u-g) \partial_n v\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. γ is the usual + * penalty parameter. */ template - void nitsche_matrix ( - FullMatrix& M, + 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; + const unsigned int n_comp = fe.get_fe().n_components(); - - Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs)); - Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs)); + 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;iF the matrices associated with the bilinear form + * @f[ + * \int_F \Bigl( \gamma [u][v] - \{\nabla u\}[v\mathbf n] - [u\mathbf + * n]\{\nabla v\} \Bigr) \; ds. + * @f] + * + * The penalty parameter should always be the mean value of the + * penalties needed for stability on each side. In the case of + * constant coefficients, it can be computed using compute_penalty(). + * + * If factor2 is missing or negative, the factor is assumed + * the same on both sides. If factors differ, note that the penalty + * parameter has to be computed accordingly. + * + * @author Guido Kanschat + * @date 2008, 2009, 2010 + */ + template + void ip_matrix ( + FullMatrix& M11, + FullMatrix& M12, + FullMatrix& M21, + FullMatrix& M22, + const FEValuesBase& fe1, + const FEValuesBase& fe2, + double penalty, + double factor1 = 1., + double factor2 = -1.) + { + const unsigned int n_dofs = fe1.dofs_per_cell; + AssertDimension(M11.n(), n_dofs); + AssertDimension(M11.m(), n_dofs); + AssertDimension(M12.n(), n_dofs); + AssertDimension(M12.m(), n_dofs); + AssertDimension(M21.n(), n_dofs); + AssertDimension(M21.m(), n_dofs); + AssertDimension(M22.n(), n_dofs); + AssertDimension(M22.m(), n_dofs); + + const double nui = factor1; + const double nue = (factor2 < 0) ? factor1 : factor2; + + for (unsigned k=0;k& n = fe1.normal_vector(k); + for (unsigned i=0;iZi the value Pi = + * pi(pi+1)/hi, where pi is + * the polynomial degree on cell Zi and + * hi is the length of Zi + * orthogonal to the current face. + * + * @author Guido Kanschat + * @date 2010 + */ + template + double compute_penalty( + const MeshWorker::DoFInfo& dinfo1, + const MeshWorker::DoFInfo& dinfo2, + unsigned int deg1, + unsigned int deg2) + { + const unsigned int normal1 = GeometryInfo::unit_normal_direction[dinfo1.face_number]; + const unsigned int normal2 = GeometryInfo::unit_normal_direction[dinfo2.face_number]; + double penalty1 = deg1 * (deg1+1) / dinfo1.cell->extent_in_direction(normal1); + double penalty2 = deg2 * (deg2+1) / dinfo2.cell->extent_in_direction(normal2); + if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children()) + { + Assert (dinfo1.face == dinfo2.face, ExcInternalError()); + Assert (dinfo1.face->has_children(), ExcInternalError()); + penalty1 *= 2; + } + const double penalty = 0.5*(penalty1 + penalty2); + return penalty; + } } } -- 2.39.5