From cf1503ac44efbb468cc173531f7637ee82c30077 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 12 Nov 2007 17:24:11 +0000 Subject: [PATCH] Implement most of the no_normal_flux function. git-svn-id: https://svn.dealii.org/trunk@15488 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-22/step-22.cc | 401 ++++++++++++++++++++++++++-- 1 file changed, 386 insertions(+), 15 deletions(-) diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index 9ddfaf893f..245db29a1f 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -43,6 +43,8 @@ #include #include #include +#include +#include #include #include @@ -139,6 +141,8 @@ double TemperatureBoundaryValues::value (const Point &p, const unsigned int /*component*/) const { +//TODO: leftover from olden times. replace by something sensible once we have +//diffusion in the temperature field if (p[0] == 0) return 1; else @@ -186,6 +190,346 @@ InitialValues::vector_value (const Point &p, +namespace +{ + /** + * A structure that stores the dim DoF + * indices that correspond to a + * vector-valued quantity at a single + * support point. + */ + template + struct VectorDoFTuple + { + unsigned int dof_indices[dim]; + + bool operator < (const VectorDoFTuple &other) const + { + for (unsigned int i=0; i other.dof_indices[i]) + return false; + return false; + } + + bool operator == (const VectorDoFTuple &other) const + { + for (unsigned int i=0; i &other) const + { + return ! (*this == other); + } + }; +} + + + +template class DH> +void +compute_no_normal_flux_constraints (const DH &dof_handler, + const unsigned int first_vector_component, + const std::set &boundary_ids, + ConstraintMatrix &constraints, + const Mapping &mapping = StaticMappingQ1::mapping) +{ + Assert (dim > 1, + ExcMessage ("This function is not useful in 1d because it amounts " + "to imposing Dirichlet values on the vector-valued " + "quantity.")); + + const FiniteElement &fe = dof_handler.get_fe(); + + std::vector face_dofs (fe.dofs_per_face); + std::vector > dof_locations (fe.dofs_per_face); + + // have a map that stores normal vectors + // for each vector-dof tuple we want to + // constrain. since we can get at the same + // vector dof tuple more than once (for + // example if it is located at a vertex + // that we visit from all adjacent cells), + // we will want to average later on the + // normal vectors computed on different + // cells as described in the documentation + // of this function. however, we can only + // average if the contributions came from + // different cells, whereas we want to + // constrain twice or more in case the + // contributions came from different faces + // of the same cell. consequently, we also + // have to store which cell a normal vector + // was computed on + typedef + std::multimap, + std::pair, typename DH::active_cell_iterator> > + DoFToNormalsMap; + + DoFToNormalsMap dof_to_normals_map; + + // now loop over all cells and all faces + typename DH::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int face_no=0; face_no < GeometryInfo::faces_per_cell; + ++face_no) + if (boundary_ids.find(cell->face(face_no)->boundary_indicator()) + != boundary_ids.end()) + { + typename DH::face_iterator face = cell->face(face_no); + + std::vector > + unit_support_points = fe.get_unit_face_support_points(); + + Quadrature aux_quad (unit_support_points); + FEFaceValues fe_values (mapping, fe, aux_quad, + update_normal_vectors); + + face->get_dof_indices (face_dofs, cell->active_fe_index()); + fe_values.reinit(cell, face_no); + + for (unsigned int i=0; i vector_dofs; + vector_dofs.dof_indices[0] = face_dofs[i]; + + for (unsigned int k=0; k= + first_vector_component) + && + (fe.face_system_to_component_index(k).first < + first_vector_component + dim)) + vector_dofs.dof_indices[fe.face_system_to_component_index(k).first] + = face_dofs[k]; + + dof_to_normals_map + .insert (std::make_pair (vector_dofs, + std::make_pair (fe_values.normal_vector(i), + cell))); + } + } + + typename DoFToNormalsMap::const_iterator + p = dof_to_normals_map.begin(); + + while (p != dof_to_normals_map.end()) + { + // first find the range of entries in + // the multimap that corresponds to the + // same vector-dof tuple. as usual, we + // define the range half-open. the + // first entry of course is 'p' + typename DoFToNormalsMap::const_iterator same_dof_range[2] + = { p }; + for (++p; p != dof_to_normals_map.end(); ++p) + if (p->first != same_dof_range[0]->first) + { + same_dof_range[1] = p; + break; + } + if (p == dof_to_normals_map.end()) + same_dof_range[1] = dof_to_normals_map.end(); + + // now compute the reverse mapping: for + // each of the cells that contributed + // to the current set of vector dofs, + // add up the normal vectors. the + // values of the map are pairs of + // normal vectors and number of cells + // that have contributed + typedef + std::map + ::active_cell_iterator, + std::pair, unsigned int> > + CellToNormalsMap; + + CellToNormalsMap cell_to_normals_map; + for (typename DoFToNormalsMap::const_iterator + q = same_dof_range[0]; + q != same_dof_range[1]; ++q) + if (cell_to_normals_map.find (q->second.second) + == cell_to_normals_map.end()) + cell_to_normals_map[q->second.second] + = std::make_pair (q->second.first, 1U); + else + { + const Tensor<1,dim> old_normal + = cell_to_normals_map[q->second.second].first; + const unsigned int old_count + = cell_to_normals_map[q->second.second].second; + + cell_to_normals_map[q->second.second] + = std::make_pair (old_normal + q->second.first, + old_count + 1); + } + + // count the maximum number of + // contributions from each cell and add + // up all the normal vectors; we will + // need the latter if each cell + // contributed exactly once -- + // otherwise we will have to compute + // something else further down below + unsigned int max_n_contributions_per_cell = 1; + Tensor<1,dim> normal; + for (typename CellToNormalsMap::const_iterator + x = cell_to_normals_map.begin(); + x != cell_to_normals_map.end(); ++x) + { + max_n_contributions_per_cell + = std::max (max_n_contributions_per_cell, + x->second.second); + normal += x->second.first; + } + + // verify that each cell can have only + // contributed at most dim times, since + // that is the maximum number of faces + // that come together at a single place + Assert (max_n_contributions_per_cell <= dim, ExcInternalError()); + + switch (max_n_contributions_per_cell) + { + // first deal with the case that a + // number of cells all have + // registered that they have a + // normal vector defined at the + // location of a given vector dof, + // and that each of them have + // encountered this vector dof + // exactly once while looping over + // all their faces. as stated in + // the documentation, this is the + // case where we want to simply + // average over all normal vectors + case 1: + { + + // compute the average normal + // vector from all the ones that + // have the same set of dofs + VectorDoFTuple dof_indices = same_dof_range[0]->first; + normal /= cell_to_normals_map.size(); + normal /= normal.norm(); + + if (cell_to_normals_map.size() == 2) + { + std::cout << "XX " << same_dof_range[0]->first.dof_indices[0] + << ' ' << same_dof_range[0]->first.dof_indices[1] + << std::endl; + std::cout << " " << cell_to_normals_map.begin()->first + << ' ' << cell_to_normals_map.begin()->second.first + << std::endl; + std::cout << " " << (++cell_to_normals_map.begin())->first + << ' ' << (++cell_to_normals_map.begin())->second.first + << std::endl; + std::cout << " " << normal << std::endl; + } + + + // then construct constraints + // from this. choose the DoF that + // has the largest component in + // the normal vector as the one + // to be constrained as this + // makes the process stable in + // cases where the normal vector + // has the form n=(1,0) or + // n=(0,1) + switch (dim) + { + case 2: + { + if (std::fabs(normal[0]) > std::fabs(normal[1])) + { + constraints.add_line (dof_indices.dof_indices[0]); + constraints.add_entry (dof_indices.dof_indices[0], + dof_indices.dof_indices[1], + -normal[1]/normal[0]); + } + else + { + constraints.add_line (dof_indices.dof_indices[1]); + constraints.add_entry (dof_indices.dof_indices[1], + dof_indices.dof_indices[0], + -normal[0]/normal[1]); + } + break; + } + + case 3: + { + if ((std::fabs(normal[0]) >= std::fabs(normal[1])) + && + (std::fabs(normal[0]) >= std::fabs(normal[2]))) + { + constraints.add_line (dof_indices.dof_indices[0]); + constraints.add_entry (dof_indices.dof_indices[0], + dof_indices.dof_indices[1], + -normal[1]/normal[0]); + constraints.add_entry (dof_indices.dof_indices[0], + dof_indices.dof_indices[2], + -normal[2]/normal[0]); + } + else + if ((std::fabs(normal[1]) >= std::fabs(normal[0])) + && + (std::fabs(normal[1]) >= std::fabs(normal[2]))) + { + constraints.add_line (dof_indices.dof_indices[1]); + constraints.add_entry (dof_indices.dof_indices[1], + dof_indices.dof_indices[0], + -normal[0]/normal[1]); + constraints.add_entry (dof_indices.dof_indices[1], + dof_indices.dof_indices[2], + -normal[2]/normal[1]); + } + else + { + constraints.add_line (dof_indices.dof_indices[2]); + constraints.add_entry (dof_indices.dof_indices[2], + dof_indices.dof_indices[0], + -normal[0]/normal[2]); + constraints.add_entry (dof_indices.dof_indices[2], + dof_indices.dof_indices[1], + -normal[1]/normal[2]); + } + + break; + } + + default: + Assert (false, ExcNotImplemented()); + } + + break; + } + + + default: + Assert (false, ExcNotImplemented()); + } + } +} + + + @@ -419,6 +763,10 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) hanging_node_constraints.clear (); DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints); + std::set no_normal_flux_boundaries; + no_normal_flux_boundaries.insert (0); + compute_no_normal_flux_constraints (dof_handler, 0, no_normal_flux_boundaries, + hanging_node_constraints); hanging_node_constraints.close (); std::vector dofs_per_component (dim+2); @@ -642,14 +990,34 @@ void BoussinesqFlowProblem::assemble_system () if (rebuild_matrices == true) { - std::vector component_mask (dim+2, true); - component_mask[dim] = component_mask[dim+1] = false; std::map boundary_values; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction(dim+2), - boundary_values, - component_mask); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + emdc = dof_handler.end(); + for (; cell!=endc; ++cell) + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (cell->vertex(v).distance(dim == 2 + ? + Point(0,-1) + : + Point(0,0,-1)) < 1e-6) + { + std::cout << "Found cell and vertex: " << cell << ' ' + << v << std::endl; + + boundary_values[cell->vertex_dof_index(v,0)] = 0; + break; + } + +// std::vector component_mask (dim+2, true); +// component_mask[dim] = component_mask[dim+1] = false; +// VectorTools::interpolate_boundary_values (dof_handler, +// 0, +// ZeroFunction(dim+2), +// boundary_values, +// component_mask); + MatrixTools::apply_boundary_values (boundary_values, system_matrix, solution, @@ -956,13 +1324,14 @@ void BoussinesqFlowProblem::assemble_rhs_T () template void BoussinesqFlowProblem::solve () { + solution = old_solution; + const InverseMatrix,SparseDirectUMFPACK> A_inverse (system_matrix.block(0,0), *A_preconditioner); Vector tmp (solution.block(0).size()); Vector schur_rhs (solution.block(1).size()); Vector tmp2 (solution.block(2).size()); - { A_inverse.vmult (tmp, system_rhs.block(0)); system_matrix.block(1,0).vmult (schur_rhs, tmp); @@ -975,15 +1344,17 @@ void BoussinesqFlowProblem::solve () SolverControl solver_control (system_matrix.block(0,0).m(), 1e-6*schur_rhs.l2_norm()); SolverCG<> cg (solver_control); - + PreconditionSSOR<> preconditioner; preconditioner.initialize (system_matrix.block(1,1), 1.2); - + + InverseMatrix,PreconditionSSOR<> > + m_inverse (system_matrix.block(1,1), preconditioner); try { cg.solve (schur_complement, solution.block(1), schur_rhs, - preconditioner); + m_inverse); } catch (...) { @@ -1172,10 +1543,10 @@ void BoussinesqFlowProblem::run () { case 2: { - GridGenerator::half_hyper_shell (triangulation, - Point(), 0.5, 1.0); + GridGenerator::hyper_shell (triangulation, + Point(), 0.5, 1.0); - static HalfHyperShellBoundary boundary; + static HyperShellBoundary boundary; triangulation.set_boundary (0, boundary); triangulation.refine_global (4); @@ -1191,7 +1562,7 @@ void BoussinesqFlowProblem::run () static HyperShellBoundary boundary; triangulation.set_boundary (0, boundary); - triangulation.refine_global (1); + triangulation.refine_global (2); break; } -- 2.39.5