-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 <int dim>
- struct VectorDoFTuple
- {
- unsigned int dof_indices[dim];
-
- bool operator < (const VectorDoFTuple<dim> &other) const
- {
- for (unsigned int i=0; i<dim; ++i)
- if (dof_indices[i] < other.dof_indices[i])
- return true;
- else
- if (dof_indices[i] > other.dof_indices[i])
- return false;
- return false;
- }
-
- bool operator == (const VectorDoFTuple<dim> &other) const
- {
- for (unsigned int i=0; i<dim; ++i)
- if (dof_indices[i] != other.dof_indices[i])
- return false;
-
- return true;
- }
-
- bool operator != (const VectorDoFTuple<dim> &other) const
- {
- return ! (*this == other);
- }
- };
- }
-}
-
-
-
-template <int dim, template <int> class DH>
-void
-compute_no_normal_flux_constraints (const DH<dim> &dof_handler,
- const unsigned int first_vector_component,
- const std::set<unsigned char> &boundary_ids,
- ConstraintMatrix &constraints,
- const Mapping<dim> &mapping = StaticMappingQ1<dim>::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<dim> &fe = dof_handler.get_fe();
-
- std::vector<unsigned int> face_dofs (fe.dofs_per_face);
- std::vector<Point<dim> > 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<dim>,
- std::pair<Tensor<1,dim>, typename DH<dim>::active_cell_iterator> >
- DoFToNormalsMap;
-
- DoFToNormalsMap dof_to_normals_map;
-
- // now loop over all cells and all faces
- typename DH<dim>::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<dim>::faces_per_cell;
- ++face_no)
- if (boundary_ids.find(cell->face(face_no)->boundary_indicator())
- != boundary_ids.end())
- {
- typename DH<dim>::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<Point<dim-1> > &
- unit_support_points = fe.get_unit_face_support_points();
- Quadrature<dim-1> aux_quad (unit_support_points);
- FEFaceValues<dim> 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<face_dofs.size(); ++i)
- if (fe.face_system_to_component_index(i).first ==
- first_vector_component)
- {
- // find corresponding other
- // components of vector
- ::internals::VectorTools::VectorDoFTuple<dim> vector_dofs;
- vector_dofs.dof_indices[0] = face_dofs[i];
-
- for (unsigned int k=0; k<fe.dofs_per_face; ++k)
- if ((k != i)
- &&
- (unit_support_points[k] == unit_support_points[i])
- &&
- (fe.face_system_to_component_index(k).first >=
- 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
- <typename DH<dim>::active_cell_iterator,
- std::pair<Tensor<1,dim>, 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<dim> &
- 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. <DIM,
- // contributions. Here we
- // do not want to average
- // the normal
- // vectors. Since we have
- // DIM contributions, let's
- // assume (and verify) that
- // they are in fact all
- // linearly independent; in
- // that case, all vector
- // components are
- // constrained and we need
- // to set them to zero
- case dim:
- {
- // assert that indeed
- // only a single cell has
- // contributed
- Assert (cell_to_normals_map.size() == 1,
- ExcInternalError());
-
- // check linear
- // independence by
- // computing the
- // determinant of the
- // matrix created from
- // all the normal
- // vectors. if they are
- // linearly independent,
- // then the determinant
- // is nonzero. if they
- // are orthogonal, then
- // the matrix is in fact
- // equal to 1 (since they
- // are all unit vectors);
- // make sure the
- // determinant is larger
- // than 1e-3 to avoid
- // cases where cells are
- // degenerate
- {
- Tensor<2,dim> t;
-
- typename DoFToNormalsMap::const_iterator x = same_dof_range[0];
- for (unsigned int i=0; i<dim; ++i, ++x)
- for (unsigned int j=0; j<dim; ++j)
- t[i][j] = x->second.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; i<dim; ++i)
- {
- constraints.add_line (same_dof_range[0]->first.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());
- }
- }
-}
-
-
-
{
cg.solve (*matrix, dst, src, preconditioner);
}
- catch (...)
+ catch (std::exception &e)
{
- Assert (false, ExcInternalError());
+ Assert (false, ExcMessage(e.what()));
}
}
hanging_node_constraints);
std::set<unsigned char> 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<unsigned int> dofs_per_component (dim+2);
{
case 2:
{
+// GridGenerator::hyper_ball (triangulation);
+
+// static const HyperBallBoundary<dim> boundary;
+// triangulation.set_boundary (0, boundary);
+
GridGenerator::hyper_cube (triangulation);
- triangulation.refine_global (5);
+ triangulation.refine_global (6);
break;
}