From 2d12212c90d823736fefbce8cb8b33cb9c8b74e1 Mon Sep 17 00:00:00 2001 From: kanschat Date: Sun, 14 Nov 2010 21:33:21 +0000 Subject: [PATCH] use LocalIntegrators instead of integrating by hand; results are still exactly the same git-svn-id: https://svn.dealii.org/trunk@22721 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-39/step-39.cc | 102 +++++++++------------------- 1 file changed, 33 insertions(+), 69 deletions(-) diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index e47b9deff3..e0cce96e4d 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -42,6 +42,11 @@ #include #include + // The include file for local + // integrators associated with the + // Laplacian +#include + // Support for multigrid methods #include #include @@ -116,47 +121,40 @@ class MatrixIntegrator : public Subscriptor // On each cell, we integrate the - // Dirichlet form. All local - // integrations consist of nested - // loops, first over all quadrature - // points and the iner loops over the - // degrees of freedom associated - // with the shape functions. + // Dirichlet form. We use the library + // of ready made integrals in + // LocalIntegrators to avoid writing + // these loops ourselves. Similarly, + // we implement Nitsche boundary + // conditions and the interior + // penalty fluxes between cells. + // + // The boundary und flux terms need a + // penalty parameter, which should be + // adjusted to the cell size and the + // polynomial degree. A safe choice + // of this parameter for constant + // coefficients can be found in + // LocalIntegrators::Laplace::compute_penalty() + // and we use this below. template void MatrixIntegrator::cell( MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationInfo& info) { - const FEValuesBase& fe = info.fe_values(); - FullMatrix& local_matrix = dinfo.matrix(0,false).matrix; - - for (unsigned int k=0; k void MatrixIntegrator::bdry( MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationInfo& info) { - const FEValuesBase& fe = info.fe_values(); - FullMatrix& local_matrix = dinfo.matrix(0,false).matrix; - - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); - - for (unsigned k=0;k::face( typename MeshWorker::IntegrationInfo& info1, typename MeshWorker::IntegrationInfo& info2) { - const FEValuesBase& fe1 = info1.fe_values(); - const FEValuesBase& fe2 = info2.fe_values(); - FullMatrix& matrix_v1u1 = dinfo1.matrix(0,false).matrix; - FullMatrix& matrix_v1u2 = dinfo1.matrix(0,true).matrix; - FullMatrix& matrix_v2u1 = dinfo2.matrix(0,true).matrix; - FullMatrix& matrix_v2u2 = dinfo2.matrix(0,false).matrix; - - const unsigned int deg = fe1.get_fe().tensor_degree(); - double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); - double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); - if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children()) - { - Assert (dinfo1.face == dinfo2.face, ExcInternalError()); - Assert (dinfo1.face->has_children(), ExcInternalError()); -// Assert (dinfo1.cell->has_children(), ExcInternalError()); - penalty1 *= 2; - } -const double penalty = penalty1 + penalty2; - - for (unsigned k=0;k