From 528c84aca088715e1edc6991c83b9b95e25aeac0 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 22 May 1998 07:40:04 +0000 Subject: [PATCH] Add projection to boundary. git-svn-id: https://svn.dealii.org/trunk@331 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/numerics/assembler.cc | 8 +- deal.II/deal.II/source/numerics/matrices.cc | 155 ++++++++++++++++++- deal.II/deal.II/source/numerics/vectors.cc | 37 ++++- 3 files changed, 190 insertions(+), 10 deletions(-) diff --git a/deal.II/deal.II/source/numerics/assembler.cc b/deal.II/deal.II/source/numerics/assembler.cc index 69c917d04e..6880086717 100644 --- a/deal.II/deal.II/source/numerics/assembler.cc +++ b/deal.II/deal.II/source/numerics/assembler.cc @@ -107,14 +107,8 @@ void Assembler::assemble (const Equation &equation) { boundary); const unsigned int n_dofs = dof_handler->get_selected_fe().total_dofs; - // clear cell matrix if (assemble_matrix) - for (unsigned int i=0; i #include #include +#include #include #include #include @@ -11,6 +12,9 @@ #include #include +#include +#include + template @@ -74,6 +78,20 @@ void MatrixCreator::create_mass_matrix (const DoFHandler &dof, +void MatrixCreator<1>::create_boundary_mass_matrix (const DoFHandler<1> &, + const FiniteElement<1> &, + const Quadrature<0> &, + const Boundary<1> &, + dSMatrix &, + const FunctionMap &, + dVector &, + vector &, + const Function<1> *) { + Assert (false, ExcNotImplemented()); +}; + + + template void MatrixCreator::create_boundary_mass_matrix (const DoFHandler &dof, const FiniteElement &fe, @@ -82,8 +100,141 @@ void MatrixCreator::create_boundary_mass_matrix (const DoFHandler & dSMatrix &matrix, const FunctionMap &rhs, dVector &rhs_vector, - vector &vec_to_dof_mapping, - const Function *a) {}; + vector &dof_to_boundary_mapping, + const Function *a) { + Assert (matrix.n() == dof.n_boundary_dofs(rhs), ExcInternalError()); + Assert (matrix.n() == matrix.m(), ExcInternalError()); + Assert (matrix.n() == rhs_vector.size(), ExcInternalError()); + Assert (rhs.size() != 0, ExcInternalError()); + Assert (dof.get_selected_fe() == fe, ExcInternalError()); + Assert (dof_to_boundary_mapping.size() == dof.n_dofs(), ExcInternalError()); + Assert (*max_element(dof_to_boundary_mapping.begin(),dof_to_boundary_mapping.end()) == + (signed int)matrix.n()-1, + ExcInternalError()); + + const unsigned int dofs_per_cell = fe.total_dofs, + dofs_per_face = fe.dofs_per_face; + dFMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + dVector cell_vector(dofs_per_cell); + + + UpdateFlags update_flags = UpdateFlags (update_JxW_values | update_q_points); + FEFaceValues fe_values (fe, q, update_flags); + + DoFHandler::active_cell_iterator cell = dof.begin_active (), + endc = dof.end (); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + // check if this face is on that part of + // the boundary we are interested in + if (rhs.find(cell->face(face)->boundary_indicator()) != rhs.end()) + { + cell_matrix.clear (); + cell_vector.clear (); + + fe_values.reinit (cell, face, fe, boundary); + + const dFMatrix &values = fe_values.get_shape_values (); + const vector &weights = fe_values.get_JxW_values (); + vector rhs_values (fe_values.n_quadrature_points); + rhs.find(cell->face(face)->boundary_indicator()) + ->second->value_list (fe_values.get_quadrature_points(), rhs_values); + + if (a != 0) + { + vector coefficient_values (fe_values.n_quadrature_points); + a->value_list (fe_values.get_quadrature_points(), coefficient_values); + for (unsigned int point=0; point dofs (dofs_per_cell); + cell->get_dof_indices (dofs); + + vector dofs_on_face_vector (dofs_per_face); + cell->face(face)->get_dof_indices (dofs_on_face_vector); + set dofs_on_face (dofs_on_face_vector.begin(), + dofs_on_face_vector.end()); + + for (unsigned int i=0; i::project_boundary_values (const DoFHandler &dof, // in this case the boundary mass matrix has // to be condensed and the solution is to // be distributed afterwards, which is not - // yet implemented + // yet implemented. The reason for this is + // that we cannot simply use the #condense# + // family of functions, since the matrices + // and vectors do not use the global + // numbering but rather the boundary + // numbering, i.e. the condense function + // needs to use another indirection. There + // should be not many technical problems, + // but it needs to be implemented if (dim<3) sparsity.compress(); else @@ -257,7 +265,34 @@ VectorTools::project_boundary_values (const DoFHandler &dof, MatrixTools::create_boundary_mass_matrix (dof, fe, q, boundary, mass_matrix, boundary_functions, rhs, dof_to_boundary_mapping); + + // same thing as above: if dim>=3 we need + // to consider constraints + Assert (dim<3, ExcNotImplemented()); + + dVector boundary_projection (rhs.size()); + + int max_iter = 1000; + double tolerance = 1.e-16; + Control control1(max_iter,tolerance); + PrimitiveVectorMemory memory(rhs.size()); + CG cg(control1,memory); + + // solve + cg (mass_matrix, boundary_projection, rhs); + + // fill in boundary values + for (unsigned int i=0; i