From: bangerth Date: Thu, 24 Jan 2008 23:28:37 +0000 (+0000) Subject: Move the stuff that computes no-flux constraints into the library proper. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8d941013f662f1147b07edec7203dc1282525159;p=dealii-svn.git Move the stuff that computes no-flux constraints into the library proper. git-svn-id: https://svn.dealii.org/trunk@15679 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc index c635db8287..fd58cac8ea 100644 --- a/deal.II/examples/step-22/step-22.cc +++ b/deal.II/examples/step-22/step-22.cc @@ -231,479 +231,6 @@ RightHandSide::vector_value (const Point &p, -namespace internals -{ - namespace VectorTools - { - /** - * 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< ::internals::VectorTools::VectorDoFTuple, - 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); - - // get the indices of the - // dofs on this cell... - face->get_dof_indices (face_dofs, cell->active_fe_index()); - - // ...and the normal - // vectors at the locations - // where they are defined: - const 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); - fe_values.reinit(cell, face_no); - - // then identify which of - // them correspond to the - // selected set of vector - // components - 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]; - - // and enter the - // (dofs,(normal_vector,cell)) - // entry into the map - dof_to_normals_map - .insert (std::make_pair (vector_dofs, - std::make_pair (fe_values.normal_vector(i), - cell))); - } - } - - // Now do something with the - // collected information. To this - // end, loop through all sets of - // pairs (dofs,normal_vector) and - // identify which entries belong to - // the same set of dofs and then do - // as described in the - // documentation, i.e. either - // average the normal vector or - // don't for this particular set of - // dofs - 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; - - Assert (old_count > 0, ExcInternalError()); - - // in the same entry, - // store again the now - // averaged normal vector - // and the new count - cell_to_normals_map[q->second.second] - = std::make_pair ((old_normal * old_count + q->second.first) / (old_count + 1), - old_count + 1); - } - - Assert (cell_to_normals_map.size() >= 1, ExcInternalError()); - - // count the maximum number of - // contributions from each cell - unsigned int max_n_contributions_per_cell = 1; - 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); - - // 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. we - // could add them up and - // divide them by the - // number of additions, - // or simply normalize - // them right away since - // we want them to have - // unit length anyway - Tensor<1,dim> normal; - for (typename CellToNormalsMap::const_iterator - x = cell_to_normals_map.begin(); - x != cell_to_normals_map.end(); ++x) - normal += x->second.first; - normal /= normal.norm(); - - const ::internals::VectorTools::VectorDoFTuple & - dof_indices = same_dof_range[0]->first; - - 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; - } - if (cell_to_normals_map.size() == 1) - { - std::cout << "YY " << 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 << " " << 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; - } - - - // this is the slightly - // more complicated case - // that a single cell has - // contributed with exactly - // DIM normal vectors to - // the same set of vector - // dofs. this is what - // happens in a corner in - // 2d and 3d (but not on an - // edge in 3d, where we - // have only 2, i.e. t; - - typename DoFToNormalsMap::const_iterator x = same_dof_range[0]; - for (unsigned int i=0; isecond.first[j]; - - Assert (std::fabs(determinant (t)) > 1e-3, - ExcMessage("Found a set of normal vectors that are nearly collinear.")); - } - - // so all components of - // this vector dof are - // constrained. enter - // this into the - // constraint matrix - for (unsigned int i=0; ifirst.dof_indices[i]); - // no add_entries here - } - - std::cout << "DIM contributions at " - << same_dof_range[0]->first.dof_indices[0] - << ' ' - << same_dof_range[0]->first.dof_indices[1] - << std::endl; - break; - } - - - // this is the case of an - // edge contribution in 3d - default: - Assert (dim >= 3, ExcInternalError()); - Assert (false, ExcNotImplemented()); - } - } -} - - - @@ -862,9 +389,9 @@ void InverseMatrix::vmult (Vector &dst, { cg.solve (*matrix, dst, src, preconditioner); } - catch (...) + catch (std::exception &e) { - Assert (false, ExcInternalError()); + Assert (false, ExcMessage(e.what())); } } @@ -939,8 +466,9 @@ void BoussinesqFlowProblem::setup_dofs (const bool setup_matrices) 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); + VectorTools::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); @@ -1727,9 +1255,14 @@ void BoussinesqFlowProblem::run () { case 2: { +// GridGenerator::hyper_ball (triangulation); + +// static const HyperBallBoundary boundary; +// triangulation.set_boundary (0, boundary); + GridGenerator::hyper_cube (triangulation); - triangulation.refine_global (5); + triangulation.refine_global (6); break; }