From: kanschat Date: Wed, 20 Oct 2010 22:20:05 +0000 (+0000) Subject: separate normal and tangential components for vector elements X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d5f968512de3804191a8590194bb03bb386359ff;p=dealii-svn.git separate normal and tangential components for vector elements git-svn-id: https://svn.dealii.org/trunk@22402 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index fb7242f7fb..25fbdcb63e 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -763,9 +763,7 @@ namespace FETools } -// This function is tested by tests/fe/internals, since it produces the matrices printed there - -//TODO:[GK] Is this correct for vector valued? +// This function is tested by tests/fe/internals, since it produces the matrices printed there template void compute_face_embedding_matrices(const FiniteElement& fe, @@ -773,29 +771,23 @@ namespace FETools const unsigned int face_coarse, const unsigned int face_fine) { + Assert(face_coarse==0, ExcNotImplemented()); + Assert(face_fine==0, ExcNotImplemented()); + const unsigned int nc = GeometryInfo::max_children_per_face; const unsigned int n = fe.dofs_per_face; const unsigned int nd = fe.n_components(); const unsigned int degree = fe.degree; - + + const bool normal = fe.conforms(FiniteElementData::Hdiv); + const bool tangential = fe.conforms(FiniteElementData::Hcurl); + for (unsigned int i=0;i tria; - GridGenerator::hyper_cube (tria, 0, 1); - tria.refine_global(1); - - MappingCartesian mapping; - QGauss q_gauss(degree+1); - const Quadrature q_fine = QProjector::project_to_face(q_gauss, face_fine); - - const unsigned int nq = q_fine.size(); - + // In order to make the loops below // simpler, we introduce vectors // containing for indices 0-n the @@ -849,6 +841,25 @@ namespace FETools } Assert (k == fe.dofs_per_face, ExcInternalError()); + // Set up meshes, one with a single + // reference cell and refine it once + Triangulation tria; + GridGenerator::hyper_cube (tria, 0, 1); + tria.refine_global(1); + MappingCartesian mapping; + + // Setup quadrature and FEValues + // for a face. We cannot use + // FEFaceValues and + // FESubfaceValues because of + // some nifty handling of + // refinement cases. Guido stops + // disliking and instead starts + // hating the anisotropic implementation + QGauss q_gauss(degree+1); + const Quadrature q_fine = QProjector::project_to_face(q_gauss, face_fine); + const unsigned int nq = q_fine.size(); + FEValues fine (mapping, fe, q_fine, update_quadrature_points | update_JxW_values | update_values); @@ -868,12 +879,21 @@ namespace FETools fine.reinit(tria.begin_active()); FullMatrix A(nq*nd, n); for (unsigned int j=0;j H(A); - + Vector v_coarse(nq*nd); Vector v_fine(n); @@ -904,10 +924,18 @@ namespace FETools // function values of the // coarse grid function in // each quadrature point. - for (unsigned int d=0;d