From 989973fdd5d700263f6f83a00192e34109cae443 Mon Sep 17 00:00:00 2001 From: guido Date: Tue, 12 Jul 2005 22:11:29 +0000 Subject: [PATCH] new function FETools::compute_face_embedding_matrices git-svn-id: https://svn.dealii.org/trunk@11135 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/fe/fe_tools.h | 35 +++++- deal.II/deal.II/source/fe/fe_tools.cc | 175 ++++++++++++++++++++++++++ 2 files changed, 206 insertions(+), 4 deletions(-) diff --git a/deal.II/deal.II/include/fe/fe_tools.h b/deal.II/deal.II/include/fe/fe_tools.h index d504f7fdd9..e79f2194df 100644 --- a/deal.II/deal.II/include/fe/fe_tools.h +++ b/deal.II/deal.II/include/fe/fe_tools.h @@ -185,11 +185,11 @@ class FETools * the finite element spaces are * actually nested. * - * @arg fe The finite element + * @param fe The finite element * class for which we compute the * embedding matrices. - * @arg matrices A pointer to - * 2dim FullMatrix + * @param matrices A pointer to + * 2dim FullMatrix * objects. This is the format * used in FiniteElementBase, * where we want to use ths @@ -198,7 +198,34 @@ class FETools template static void compute_embedding_matrices(const FiniteElement &fe, FullMatrix* matrices); - + + /** + * Compute the embedding matrices + * on faces needed for constraint + * matrices. + * + * @param fe The finite element + * for which to compute these + * matrices. + * @param matrices An array of + * 2dim-1 FullMatrix + * objects,holding the embedding + * matrix for each subface. + * @param face_coarse The number + * of the face on the coarse side + * of the face for which this is + * computed. + * @param face_fine The number + * of the face on the refined side + * of the face for which this is + * computed. + */ + template + static void compute_face_embedding_matrices(const FiniteElement& fe, + FullMatrix* matrices, + unsigned int face_coarse, + unsigned int face_fine); + /** * Compute the * L2-projection diff --git a/deal.II/deal.II/source/fe/fe_tools.cc b/deal.II/deal.II/source/fe/fe_tools.cc index 41e682a1cd..e6285240b3 100644 --- a/deal.II/deal.II/source/fe/fe_tools.cc +++ b/deal.II/deal.II/source/fe/fe_tools.cc @@ -13,6 +13,7 @@ #include +#include #include #include #include @@ -586,6 +587,175 @@ FETools::compute_embedding_matrices(const FiniteElement& fe, } +// This function is tested by tests/fe/internals, since it produces the matrices printed there + +//TODO:[GK] Is this correct for vector valued? +template +void +FETools::compute_face_embedding_matrices(const FiniteElement& fe, + FullMatrix* matrices, + unsigned int face_coarse, + unsigned int face_fine) +{ + const unsigned int nc = GeometryInfo::subfaces_per_face; + const unsigned int n = fe.dofs_per_face; + const unsigned int nd = fe.n_components(); + const unsigned int degree = fe.degree; + + 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.n_quadrature_points; + + // In order to make the loops below + // simpler, we introduce vectors + // containing for indices 0-n the + // number of the corresponding + // shape value on the cell. + std::vector face_c_dofs(n); + std::vector face_f_dofs(n); + unsigned int k=0; + for (unsigned int i=0;i::vertices_per_face;++i) + { + const unsigned int offset_c = GeometryInfo::face_to_cell_vertices(face_coarse, i) + *fe.dofs_per_vertex; + const unsigned int offset_f = GeometryInfo::face_to_cell_vertices(face_fine, i) + *fe.dofs_per_vertex; + for (unsigned int j=0;j::lines_per_face;++i) + { + const unsigned int offset_c = fe.first_line_index + + GeometryInfo::face_to_cell_lines(face_coarse, i-1) + *fe.dofs_per_line; + const unsigned int offset_f = fe.first_line_index + + GeometryInfo::face_to_cell_lines(face_fine, i-1) + *fe.dofs_per_line; + for (unsigned int j=0;j::quads_per_face;++i) + { + const unsigned int offset_c = fe.first_quad_index + + face_coarse + *fe.dofs_per_quad; + const unsigned int offset_f = fe.first_quad_index + + face_fine + *fe.dofs_per_quad; + for (unsigned int j=0;j fine (mapping, fe, q_fine, + update_q_points | update_JxW_values | update_values); + + // We search for the polynomial on + // the small cell, being equal to + // the coarse polynomial in all + // quadrature points. + + // First build the matrix for this + // least squares problem. This + // contains the values of the fine + // cell polynomials in the fine + // cell grid points. + + // This matrix is the same for all + // children. + fine.reinit(tria.begin_active()); + FullMatrix A(nq*nd, n); + for (unsigned int d=0;d H(A); + + Vector v_coarse(nq*nd); + Vector v_fine(n); + + + + for (unsigned int cell_number = 0; cell_number < GeometryInfo::subfaces_per_face; + ++cell_number) + { + const Quadrature q_coarse + = QProjector::project_to_subface(q_gauss, face_coarse, cell_number); + FEValues coarse (mapping, fe, q_coarse, update_values); + + typename Triangulation::active_cell_iterator fine_cell + = tria.begin(0)->child(GeometryInfo::child_cell_on_face(face_coarse, + cell_number)); + fine.reinit(fine_cell); + coarse.reinit(tria.begin(0)); + + FullMatrix &this_matrix = matrices[cell_number]; + + // Compute this once for each + // coarse grid basis function + for (unsigned int i=0;i void @@ -1553,6 +1723,11 @@ template void FETools::compute_embedding_matrices (const FiniteElement &, FullMatrix*); +template +void FETools::compute_face_embedding_matrices +(const FiniteElement &, FullMatrix*, + unsigned int, unsigned int); + template void FETools::compute_projection_matrices (const FiniteElement &, FullMatrix*); -- 2.39.5