From: Guido Kanschat Date: Thu, 5 Jan 2017 17:07:23 +0000 (+0100) Subject: tangential jumps for Laplacian X-Git-Tag: v8.5.0-rc1~276^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3b48b751ee4a72edc0ac166b6861f9f750a878bd;p=dealii.git tangential jumps for Laplacian --- diff --git a/include/deal.II/integrators/laplace.h b/include/deal.II/integrators/laplace.h index 7379965b13..e7a156eac1 100644 --- a/include/deal.II/integrators/laplace.h +++ b/include/deal.II/integrators/laplace.h @@ -184,6 +184,65 @@ namespace LocalIntegrators } } + /** + * Weak boundary condition of Nitsche type for the Laplacian applied to the + * tangential component only, namely on + * the face F the matrix + * @f[ + * \int_F \Bigl(\gamma u_\tau v_\tau - \partial_n u_\tau v_\tau - u_\tau \partial_n v_\tau\Bigr)\;ds. + * @f] + * + * Here, $\gamma$ is the penalty parameter suitably computed with + * compute_penalty(). + * + * @author Guido Kanschat + * @date 2017 + */ + template + void nitsche_tangential_matrix ( + FullMatrix &M, + const FEValuesBase &fe, + double penalty, + double factor = 1.) + { + const unsigned int n_dofs = fe.dofs_per_cell; + AssertDimension(fe.get_fe().n_components(), dim); + Assert (M.m() == n_dofs, ExcDimensionMismatch(M.m(), n_dofs)); + Assert (M.n() == n_dofs, ExcDimensionMismatch(M.n(), n_dofs)); + + for (unsigned int k=0; k n = fe.normal_vector(k); + for (unsigned int i=0; iF the vector @@ -428,7 +487,7 @@ namespace LocalIntegrators ngradv2n += n*fe2.shape_grad_component(i,k,d)*n[d]; } - for (unsigned int d=0; d +#include +#include + +#include +#include +#include +#include + +using namespace LocalIntegrators::Laplace; + + + +template +void test_boundary(const FiniteElement &fe, bool diff=false) +{ + Triangulation tr; + GridGenerator::hyper_cube(tr); + + DoFHandler dof(tr); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + DoFTools::make_zero_boundary_constraints(dof, constraints); + + QGauss quadrature(fe.tensor_degree()+1); + FEFaceValues fev(fe, quadrature, update_values | update_gradients | update_normal_vectors | update_JxW_values); + + FullMatrix M(fe.dofs_per_cell); + FullMatrix Mglobal(dof.n_dofs()); + std::vector indices(fe.dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = dof.begin_active(); + cell->get_dof_indices(indices); + for (unsigned i=0; i::faces_per_cell; ++i) + { + fev.reinit(cell, i); + nitsche_tangential_matrix(M, fev, 10.); + if (diff) + nitsche_matrix(M, fev, 10., -1.); + constraints.distribute_local_to_global(M, indices, indices, Mglobal); + } + deallog << fe.get_name() << ": bdry norm " << Mglobal.frobenius_norm() << std::endl; +} + + + +template +void test_face(const FiniteElement &fe, bool diff=false) +{ + Triangulation tr; + GridGenerator::hyper_cube(tr); + tr.refine_global(); + + DoFHandler dof(tr); + dof.distribute_dofs(fe); + ConstraintMatrix constraints; + DoFTools::make_zero_boundary_constraints(dof, constraints); + + QGauss quadrature(fe.tensor_degree()+1); + FEFaceValues fev1(fe, quadrature, update_values | update_gradients | update_normal_vectors | update_JxW_values); + FEFaceValues fev2(fe, quadrature, update_values | update_gradients | update_normal_vectors | update_JxW_values); + + FullMatrix M11(fe.dofs_per_cell); + FullMatrix M12(fe.dofs_per_cell); + FullMatrix M21(fe.dofs_per_cell); + FullMatrix M22(fe.dofs_per_cell); + FullMatrix Mglobal(dof.n_dofs()); + std::vector indices1(fe.dofs_per_cell); + std::vector indices2(fe.dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell1 = dof.begin_active(); + typename DoFHandler::active_cell_iterator cell2 = ++dof.begin_active(); + + cell1->get_dof_indices(indices1); + cell2->get_dof_indices(indices2); + + fev1.reinit(cell1, 1); + fev2.reinit(cell2, 0); + ip_tangential_matrix(M11, M12, M21, M22, fev1, fev2, 10); + + if (diff) + ip_matrix(M11, M12, M21, M22, fev1, fev2, 10, -1.); + + constraints.distribute_local_to_global(M11, indices1, indices1, Mglobal); + constraints.distribute_local_to_global(M21, indices2, indices1, Mglobal); + constraints.distribute_local_to_global(M12, indices1, indices2, Mglobal); + constraints.distribute_local_to_global(M22, indices2, indices2, Mglobal); + deallog << fe.get_name() << ": face norm " << Mglobal.frobenius_norm() << std::endl; +} + + + +template +void +test() +{ + // In each dimension, the first four outputs should be zero: the + // tangential jump of Nedelec elements is zero and for + // Raviart-Thomas, the tangential and the full jump are the same. + FE_Nedelec n1(2); + test_boundary(n1); + test_face(n1); + + FE_RaviartThomas r1(2); + test_boundary(r1, true); + test_face(r1,true); + test_boundary(r1); + test_face(r1); + + FE_DGQ q1(1); + FESystem sys1(q1,dim); + test_boundary(sys1); + test_boundary(sys1, true); + test_face(sys1); + test_face(sys1, true); + deallog << std::endl; +} + + +int main() +{ + initlog(); + deallog.threshold_double(1.e-10); + + test<2>(); + test<3>(); +} diff --git a/tests/integrators/laplacian_02.output b/tests/integrators/laplacian_02.output new file mode 100644 index 0000000000..0024ef3b7f --- /dev/null +++ b/tests/integrators/laplacian_02.output @@ -0,0 +1,23 @@ + +DEAL::FE_Nedelec<2>(2): bdry norm 0 +DEAL::FE_Nedelec<2>(2): face norm 0 +DEAL::FE_RaviartThomas<2>(2): bdry norm 0 +DEAL::FE_RaviartThomas<2>(2): face norm 0 +DEAL::FE_RaviartThomas<2>(2): bdry norm 1143.78 +DEAL::FE_RaviartThomas<2>(2): face norm 259.355 +DEAL::FESystem<2>[FE_DGQ<2>(1)^2]: bdry norm 52.2707 +DEAL::FESystem<2>[FE_DGQ<2>(1)^2]: bdry norm 52.2707 +DEAL::FESystem<2>[FE_DGQ<2>(1)^2]: face norm 4.28174 +DEAL::FESystem<2>[FE_DGQ<2>(1)^2]: face norm 4.28174 +DEAL:: +DEAL::FE_Nedelec<3>(2): bdry norm 0 +DEAL::FE_Nedelec<3>(2): face norm 0 +DEAL::FE_RaviartThomas<3>(2): bdry norm 0 +DEAL::FE_RaviartThomas<3>(2): face norm 0 +DEAL::FE_RaviartThomas<3>(2): bdry norm 5869.16 +DEAL::FE_RaviartThomas<3>(2): face norm 1270.58 +DEAL::FESystem<3>[FE_DGQ<3>(1)^3]: bdry norm 86.7436 +DEAL::FESystem<3>[FE_DGQ<3>(1)^3]: bdry norm 47.9857 +DEAL::FESystem<3>[FE_DGQ<3>(1)^3]: face norm 1.59571 +DEAL::FESystem<3>[FE_DGQ<3>(1)^3]: face norm 1.12834 +DEAL::