* under consideration has full dimension, i.e. constitutes a cell.
*
* @ingroup dofs
- * @ingroup Accessors
+ * @ingroup Accessors
* @author Wolfgang Bangerth, 1999
*/
template <int structdim, int dim, int spacedim>
- struct Inheritance
+ struct Inheritance
{
/**
* Declaration of the @p typedef.
* since they are more secure to changes in the class naming and
* template interface as well as providing easier typing (much less
* complicated names!).
- *
+ *
*
* <h3>Inheritance</h3>
*
/**
* @{
*/
-
+
/**
* Default constructor. Provides
* an accessor that can't be
* used.
*/
DoFAccessor ();
-
+
/**
* Constructor
*/
/**
* @}
*/
-
+
/**
* Return a handle on the
* DoFHandler object which we
*/
DoFAccessor<structdim,DH> &
operator = (const DoFAccessor<structdim,DH> &da);
-
+
/**
* Implement the copy operator needed
* for the iterator classes.
* coordinates of the TriaAccessor.
*/
void copy_from (const TriaAccessorBase<structdim, DH::dimension, DH::space_dimension> &da);
-
+
/**
* @name Accessing sub-objects
*/
/**
* @{
*/
-
+
/**
* Return the indices of the dofs of this
* object in the standard ordering: dofs
/**
* @{
*/
-
+
/**
* Return the number of finite
* elements that are active on a
*/
unsigned int
nth_active_fe_index (const unsigned int n) const;
-
+
/**
* Return true if the finite
* element with given index is
/**
* @}
*/
-
+
/**
* Exceptions for child classes
*
DeclException0 (ExcNotActive);
/**
* Exception
- *
+ *
* @ingroup Exceptions
*/
DeclException0 (ExcCantCompareIterators);
- protected:
+ protected:
/**
* Store the address of the DoFHandler object
DH *dof_handler;
/**
- * Compare for equality.
+ * Compare for equality.
*/
bool operator == (const DoFAccessor &) const;
-
+
/**
* Compare for inequality.
*/
const unsigned int i,
const unsigned int index,
const unsigned int fe_index = DH::default_fe_index) const;
-
+
/**
* Iterator classes need to be friends
* because they need to access operator==
* Extract space dimension from DH.
*/
static const unsigned int spacedim = DH::space_dimension;
-
+
/**
* Declare the data type that
* this accessor class expects to
/**
* @{
*/
-
+
/**
* Return the @p ith neighbor as
* a DoF cell iterator. This
* of this cell.
*
* This function is not implemented in
- * 1D, and maps to DoFAccessor<2,
- * dim>::line in 2D.
+ * 1D, and maps to DoFAccessor::line
+ * in 2D.
*/
typename internal::DoFHandler::Iterators<DH>::face_iterator
face (const unsigned int i) const;
template <class OutputVector, typename number>
void set_dof_values (const Vector<number> &local_values,
OutputVector &values) const;
-
+
/**
* Return the interpolation of
* the given finite element
void
distribute_local_to_global (const FullMatrix<number> &local_source,
OutputMatrix &global_destination) const;
-
+
/**
* This function does what the two
* <tt>distribute_local_to_global</tt>
const Vector<number> &local_vector,
OutputMatrix &global_matrix,
OutputVector &global_vector) const;
-
+
/**
* @}
*/
/**
* @{
*/
-
+
/**
* Return the indices of the dofs of this
* quad in the standard ordering: dofs
/**
* @}
*/
-
+
private:
/**
* Update the cache in which we
// $Id$
// Version: $Name$
//
-// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+// Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* version 5.2. That is, the vertex and face ordering in 2d is assumed
* to be
* @verbatim
- * 2
+ * 2
* 3--->---2
* | |
* 3^ ^1
* Triangulation::create_triangulation_compatibility() (and not the
* Triangulation::create_triangulation()) function must be used. For
* a typical use of the reorder_cells() function see the
- * implementation of the GridIn <tt>::read_*()</tt> functions.
+ * implementation of the GridIn <code>read_*()</code> functions.
*
*
* <h3>Statement of problems</h3>
* @verbatim
* o
* / \
- * o o
+ * o o
* | \ / |
- * | o |
+ * | o |
* | | |
* o---o---o
* @endverbatim
* directions; it will soon be obvious that there can't exists such a
* thing, even if we allow that there might be cells with clockwise
* and counterclockwise orientation of the lines at the same time.)
- *
+ *
* One might argue that the definition of unique directions for faces
* and edges, and the definition of directions relative to the cells
* they bound, is a misfeature of deal.II. In fact, it makes reading
* ends of the string of cells, the string is twisted by 180 degrees,
* then no such orientation is possible any more, as can easily be
* checked. In effect, some meshes could not be used in deal.II.
- * In order to overcome this problem, the <code>face_rotation</code>,
+ * In order to overcome this problem, the <code>face_rotation</code>,
* <code>face_flip</code> and <code>line_orientation</code> flags have
* been introduced. With these, it is possible to treat all purely hexahedral
* meshes. However, in order to reduce the effect of possible bugs, it should
* 0->-1------2
* @endverbatim
*
- * The point here is the following: assume we wanted to prolong the grid to
+ * The point here is the following: assume we wanted to prolong the grid to
* the left like this:
* @verbatim
* o---o---o---o---o------o
* implemented in 2000 by Wolfgang Bangerth. Although it is no longer used
* here is how it works, and why it doesn't always work for large meshes since
* its run-time can exponential in bad cases.
- *
+ *
* The first observation is that although there are counterexamples,
* problems are usually local. For example, in the second example
* mentioned above, if we had numbered the cells in a way that
class GridReordering
{
public:
-
+
/**
* This is the main function,
* doing what is announced in
*/
static void invert_all_cells_of_negative_grid(
const std::vector<Point<spacedim> > &all_vertices,
- std::vector<CellData<dim> > &original_cells);
+ std::vector<CellData<dim> > &original_cells);
};
* @code
* i->set_refine_flag ();
* @endcode
- *
+ *
* Iterators are used whenever a loop over all lines, quads, cells
* etc. is to be performed. These loops can then be coded like this:
* @code
* a fixed value <tt>end</tt> inside the loop instead of
* <tt>tria.end()</tt>, since the creation and copying of these
* iterators is rather expensive compared to normal pointers.
- *
+ *
* The objects pointed to are accessors, derived from
* TriaAccessorBase. Which kind of accessor is determined by the template
* argument <em>Accessor</em>. These accessors are not so much data
* functions from the algorithm section of the C++ standard,
* e.g. <em>count_if</em> (see the documentation for Triangulation for
* an example) and several others.
- *
+ *
* <h3>Implementation</h3>
*
* The iterator class itself does not have much functionality. It only
* TriaAccessor, DoFAccessor and MGDoFAccessor,
* respectively. In each group, there is an accessor to cells, which
* have more functionality.
- *
+ *
* @attention It seems impossible to preserve constness of a
* triangulation through iterator usage. Thus, if you declare pointers
* to a <tt>const</tt> triangulation object, you should be well aware
* in the documentation of TriaAccessorBase, where the iterator states are
* checked and implemented.
*
- *
+ *
* <h3>Past-the-end iterators</h3>
- *
+ *
* There is a representation of past-the-end-pointers, denoted by special
* values of the member variables @p present_level and @p present_index:
* If <tt>present_level>=0</tt> and <tt>present_index>=0</tt>, then the object is valid
* Past-the-end iterators may also be used to compare an iterator with the
* <i>before-the-start</i> value, when running backwards. There is no
* distiction between the iterators pointing past the two ends of a vector.
- *
+ *
* By defining only one value to be past-the-end and making all other values
* invalid provides a second track of security: if we should have forgotten
* a check in the library when an iterator is incremented or decremented,
*/
template <typename Accessor>
class TriaRawIterator :
-#ifdef HAVE_STD_ITERATOR_CLASS
+#ifdef HAVE_STD_ITERATOR_CLASS
public std::iterator<std::bidirectional_iterator_tag,Accessor>
#else
public bidirectional_iterator<Accessor,int>
* implementation actually is.
*/
typedef Accessor AccessorType;
-
+
/**
* Empty constructor. Such an object
* is not usable!
*/
template <typename OtherAccessor>
TriaRawIterator (const OtherAccessor &a);
-
+
/**
* Proper constructor, initialized
* with the triangulation, the
* class of this object. One such path
* would be derived class to base class,
* which for example may be used to get
- * a Triangulation@p ::raw_cell_iterator from
- * a DoFHandler@p ::raw_cell_iterator, since
+ * a Triangulation::raw_cell_iterator from
+ * a DoFHandler::raw_cell_iterator, since
* the DoFAccessor class is derived from
* the TriaAccessorBase class.
*/
*/
template <typename OtherAccessor>
TriaRawIterator (const TriaActiveIterator<OtherAccessor> &i);
-
+
/**
* @name Dereferencing
*/
* past the end iterators.
*/
const Accessor & operator * () const;
-
+
/**
* Dereferencing operator, non-@p const
* version.
*/
Accessor & operator * ();
-
+
/**
* Dereferencing operator, returns a
* reference of the cell pointed to.
* version.
*/
const Accessor * operator -> () const;
-
+
/**
* Dereferencing operator, non-@p const
* version.
*/
Accessor * operator -> ();
/*@}*/
-
+
/**
* Assignment operator.
*/
TriaRawIterator & operator = (const TriaRawIterator &);
-
+
/**
* Assignment operator.
*/
// template <class OtherAccessor>
// TriaRawIterator & operator = (const TriaRawIterator<OtherAccessor>&);
-
+
/**
* Assignment operator.
*/
// template <class OtherAccessor>
// TriaRawIterator & operator = (const TriaIterator<OtherAccessor>&);
-
+
/**
* Assignment operator.
*/
// template <class OtherAccessor>
// TriaRawIterator & operator = (const TriaActiveIterator<OtherAccessor>&);
-
+
/**
* Compare for equality.
*/
bool operator == (const TriaRawIterator &) const;
-
+
/**
* Compare for inequality.
*/
* is returned in that case.
*/
bool operator < (const TriaRawIterator &) const;
-
+
/**@name Advancement of iterators*/
/*@{*/
/**
* faces have no level.
*/
TriaRawIterator & operator ++ ();
-
+
/**
* Postfix <tt>++</tt> operator: <tt>i++</tt>. This
* operator advances the iterator to
* faces have no level.
*/
TriaRawIterator & operator -- ();
-
+
/**
* Postfix @p -- operator: @p i--. This
* operator advances the iterator to
*/
TriaRawIterator operator -- (int);
/*@}*/
-
+
/**
* Return the state of the iterator.
*/
<< (arg1.state() == IteratorState::valid ? "valid" :
(arg1.state() == IteratorState::past_the_end ?
"past_the_end" : "invalid")));
-
+
/**
* Exception
*/
* Exception
*/
DeclException0 (ExcInvalidComparison);
-
+
/*@}*/
protected:
/**
* @ingroup Iterators
*/
template <typename Accessor>
-class TriaIterator : public TriaRawIterator<Accessor>
+class TriaIterator : public TriaRawIterator<Accessor>
{
public:
/**
const int level,
const int index,
const typename Accessor::AccessorData *local_data = 0);
-
+
/**
* This is a conversion operator
* (constructor) which takes another
* class of this object. One such path
* would be derived class to base class,
* which for example may be used to get
- * a Triangulation@p ::cell_iterator from
- * a DoFHandler@p ::cell_iterator, since
+ * a Triangulation::cell_iterator from
+ * a DoFHandler::cell_iterator, since
* the DoFAccessor class is derived from
* the TriaAccessorBase class.
*/
*/
template <typename OtherAccessor>
TriaIterator (const TriaActiveIterator<OtherAccessor> &i);
-
+
/**
* Assignment operator.
*/
TriaIterator<Accessor> &
operator = (const TriaIterator<Accessor>&);
-
+
/**
* Cross assignment operator. This
* assignment is only valid if the
*/
TriaIterator<Accessor> &
operator = (const TriaRawIterator<Accessor>&);
-
+
/**
* Assignment
* operator. Requires, that
*/
TriaIterator<Accessor> operator -- (int);
/*@}*/
-
+
/**
* Exception
*/
* @ingroup Iterators
*/
template <typename Accessor>
-class TriaActiveIterator : public TriaIterator<Accessor>
+class TriaActiveIterator : public TriaIterator<Accessor>
{
public:
/**
* class of this object. One such path
* would be derived class to base class,
* which for example may be used to get
- * a Triangulation@p ::active_cell_iterator from
- * a DoFHandler@p ::active_cell_iterator, since
+ * a Triangulation::active_cell_iterator from
+ * a DoFHandler::active_cell_iterator, since
* the DoFAccessor class is derived from
* the TriaAccessorBase class.
*/
*/
TriaActiveIterator<Accessor> operator -- (int);
/*@}*/
-
+
/**
* Exception
*/
Assert (Accessor::structure_dimension==Accessor::dimension ||
state() == IteratorState::valid,
ExcDereferenceInvalidObject(accessor));
-
+
return accessor;
}
// has_children() is called anyway, even if
// state==IteratorState::past_the_end, and will then
// throw the exception!
- if (this->state() != IteratorState::past_the_end)
+ if (this->state() != IteratorState::past_the_end)
Assert (this->accessor.has_children()==false,
ExcAssignmentOfInactiveObject());
#endif
// $Id$
// Version: $Name$
//
-// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* want to add several matrices and could then condense afterwards
* only once, instead of for every matrix. To actually do computations
* with these matrices, you have to condense the matrix using the
- * ConstraintMatrix@p ::condense function; you also have to
+ * ConstraintMatrix::condense function; you also have to
* condense the right hand side accordingly and distribute the
* solution afterwards.
*
* top of a corresponding sparsity pattern. Since we only work on a subset of
* the degrees of freedom, we can't use the matrices and sparsity patterns
* that are created for the entire set of degrees of freedom. Rather, you
- * should use the DoFHandler@p ::make_boundary_sparsity_pattern() function to
+ * should use the DoFHandler::make_boundary_sparsity_pattern() function to
* create the correct sparsity pattern, and build a matrix on top of it.
- *
+ *
* Note that at present there is no function that computes the mass matrix for
* <em>all</em> shape functions, though such a function would be trivial to
* implement.
{
public:
/**
- * Assemble the mass matrix. If no
+ * Assemble the mass matrix. If no
* coefficient is given, it is assumed
* to be unity.
*
Vector<double> &rhs_vector,
const Function<spacedim> * const a = 0);
-
+
/**
* Assemble the mass matrix and a
* right hand side vector along
* matrix. If no coefficient is
* given, it is assumed to be
* constant one.
- *
+ *
* If the library is configured
* to use multithreading, this
* function works in parallel.
const Quadrature<dim> &q,
SparseMatrix<double> &matrix,
const Function<spacedim> * const a = 0);
-
+
/**
* Calls the
* create_laplace_matrix()
* vector. If no coefficient is
* given, it is assumed to be
* constant one.
- *
+ *
* If the library is configured
* to use multithreading, this
* function works in parallel.
const hp::QCollection<dim> &q,
SparseMatrix<double> &matrix,
const Function<spacedim> * const a = 0);
-
+
/**
* Like the functions above, but for hp
* dof handlers, mappings, and quadrature
/**
* Exception
*/
- DeclException0 (ExcComponentMismatch);
+ DeclException0 (ExcComponentMismatch);
private:
/**
* Convenience abbreviation for
* dof handler that shouls be used.
*/
template <typename DH>
- struct IteratorRange
+ struct IteratorRange
{
/**
* Typedef for the iterator type.
* iterators.
*/
typedef std::pair<active_cell_iterator,active_cell_iterator> iterator_pair;
-
+
/**
* Constructor. Initialize
* the two values by the
* initialization.
*/
IteratorRange (const iterator_pair &ip);
-
+
/**
* Pair of iterators denoting
* a half-open range.
*/
active_cell_iterator first, second;
};
-
+
/**
* Version of the same function
* (without suffix @p _1) with
*
* There are two ways to incorporate fixed degrees of freedom such as boundary
* nodes into a linear system, as discussed below.
- *
+ *
*
* <h3>Global elimination</h3>
*
* processed node is zero, the Gauss step does not change the right
* hand side. We need therefore not take special care of other
* boundary nodes.
- *
+ *
* To make solving faster, we preset the solution vector with the right
* boundary values (as to why this is necessary, see the discussion below in
* the description of local elimination). It it not clear whether the deletion
* In conclusion, local elimination of boundary nodes only works if
* there are no hanging nodes and even then doesn't always work fully
* satisfactorily.
- *
+ *
* @ingroup numerics
* @author Wolfgang Bangerth, 1998, 2000, 2004, 2005
*/
* documentation. This function works on
* the classes that are used to wrap
* PETSc objects.
- *
+ *
* Note that this function is not very
* efficient: it needs to alternatingly
* read and write into the matrix, a
PETScWrappers::Vector &solution,
PETScWrappers::MPI::Vector &right_hand_side,
const bool eliminate_columns = true);
-
+
/**
* Same as above but for BlockSparseMatrix.
- */
+ */
static void
apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
PETScWrappers::MPI::BlockSparseMatrix &matrix,
PETScWrappers::MPI::BlockVector &solution,
PETScWrappers::MPI::BlockVector &right_hand_side,
const bool eliminate_columns = true);
-
+
#endif
#ifdef DEAL_II_USE_TRILINOS
* function works on the classes
* that are used to wrap Trilinos
* objects.
- *
+ *
* Note that this function is not
* very efficient: it needs to
* alternatingly read and write
TrilinosWrappers::Vector &solution,
TrilinosWrappers::Vector &right_hand_side,
const bool eliminate_columns = true);
-
+
/**
* This function does the same as
* the one above, except now
* function works on the classes
* that are used to wrap Trilinos
* objects.
- *
+ *
* Note that this function is not
* very efficient: it needs to
* alternatingly read and write
TrilinosWrappers::MPI::Vector &solution,
TrilinosWrappers::MPI::Vector &right_hand_side,
const bool eliminate_columns = true);
-
+
/**
* This function does the same as
* the one above, except now working
TrilinosWrappers::MPI::BlockVector &right_hand_side,
const bool eliminate_columns = true);
#endif
-
+
/**
* Rather than applying boundary
* values to the global matrix
FullMatrix<double> &local_matrix,
Vector<double> &local_rhs,
const bool eliminate_columns);
-
+
/**
* Exception
*/
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* vector is reinitialized to the new size and filled with the interpolated
* values. This class is used in the step-15,
* step-31, and step-33 tutorial programs.
- *
+ *
* <h3>Usage</h3>
*
* As the interpolation while
- * coarsening is much more complicated to organize
+ * coarsening is much more complicated to organize
* (see further documentation below) than interpolation while pure refinement,
* @p SolutionTransfer offers two possible usages.
* <ul>
* soltrans.interpolate(solution, interpolated_solution);
* @endverbatim
*
- * Multiple calls to the function
+ * Multiple calls to the function
* <tt>interpolate (const Vector<number> &in, Vector<number> &out)</tt>
* are NOT allowed. Interpolating several functions can be performed in one step
- * by using
+ * by using
* <tt>void interpolate (const vector<Vector<number> >&all_in, vector<Vector<number> >&all_out) const</tt>,
- * and using the respective @p prepare_for_coarsening_and_refinement function
+ * and using the respective @p prepare_for_coarsening_and_refinement function
* taking several vectors as input before actually refining and coarsening the
* triangulation (see there).
* </ul>
* discrete function @p in interpolated on the refined mesh.
*
* The <tt>refine_interpolate(in,out)</tt> function can be called multiple times for
- * arbitrary many discrete functions (solution vectors) on the original grid.
+ * arbitrary many discrete functions (solution vectors) on the original grid.
*
- * <li> Solution transfer while coarsening and refinement. After
- * calling Triangulation@p ::prepare_coarsening_and_refinement the
- * coarsen flags of either all or none of the children of a
- * (father-)cell are set. While coarsening
- * (Triangulation@p ::execute_coarsening_and_refinement)
+ * <li> Solution transfer while coarsening and refinement. After
+ * calling Triangulation::prepare_coarsening_and_refinement the
+ * coarsen flags of either all or none of the children of a
+ * (father-)cell are set. While coarsening
+ * (Triangulation::execute_coarsening_and_refinement)
* the cells that are not needed any more will be deleted from the Triangulation.
- *
+ *
* For the interpolation from the (to be coarsenend) children to their father
* the children cells are needed. Hence this interpolation
* and the storing of the interpolated values of each of the discrete functions
* that we want to interpolate needs to take place before these children cells
* are coarsened (and deleted!!). Again a pointers for the relevant cells is
- * set to point to these values (see below).
+ * set to point to these values (see below).
* Additionally the DoF indices of the cells
* that will not be coarsened need to be stored according to the solution
* transfer while pure refinement (cf there). All this is performed by
- * <tt>prepare_for_coarsening_and_refinement(all_in)</tt> where the
+ * <tt>prepare_for_coarsening_and_refinement(all_in)</tt> where the
* <tt>vector<Vector<number> >vector all_in</tt> includes
* all discrete functions to be interpolated onto the new grid.
*
* As we need two different kinds of pointers (<tt>vector<unsigned int> *</tt> for the Dof
* indices and <tt>vector<Vector<number> > *</tt> for the interpolated DoF values)
* we use the @p Pointerstruct that includes both of these pointers and
- * the pointer for each cell points to these @p Pointerstructs.
- * On each cell only one of the two different pointers is used at one time
+ * the pointer for each cell points to these @p Pointerstructs.
+ * On each cell only one of the two different pointers is used at one time
* hence we could use a
- * <tt>void * pointer</tt> as <tt>vector<unsigned int> *</tt> at one time and as
+ * <tt>void * pointer</tt> as <tt>vector<unsigned int> *</tt> at one time and as
* <tt>vector<Vector<number> > *</tt> at the other but using this @p Pointerstruct
* in between makes the use of these pointers more safe and gives better
* possibility to expand their usage.
- *
+ *
* In <tt>interpolate(all_in, all_out)</tt> the refined cells are treated according
* to the solution transfer while pure refinement. Additionally, on each
- * cell that is coarsened (hence previously was a father cell),
+ * cell that is coarsened (hence previously was a father cell),
* the values of the discrete
- * functions in @p all_out are set to the stored local interpolated values
- * that are accessible due to the 'vector<Vector<number> > *' pointer in
+ * functions in @p all_out are set to the stored local interpolated values
+ * that are accessible due to the 'vector<Vector<number> > *' pointer in
* @p Pointerstruct that is pointed to by the pointer of that cell.
* It is clear that <tt>interpolate(all_in, all_out)</tt> only can be called with
* the <tt>vector<Vector<number> > all_in</tt> that previously was the parameter
- * of the <tt>prepare_for_coarsening_and_refinement(all_in)</tt> function. Hence
- * <tt>interpolate(all_in, all_out)</tt> can (in contrast to
+ * of the <tt>prepare_for_coarsening_and_refinement(all_in)</tt> function. Hence
+ * <tt>interpolate(all_in, all_out)</tt> can (in contrast to
* <tt>refine_interpolate(in, out)</tt>) only be called once.
* </ul>
*
class SolutionTransfer
{
public:
-
+
/**
* Constructor, takes the current DoFHandler
* as argument.
* Destructor
*/
~SolutionTransfer();
-
+
/**
* Reinit this class to the state that
* it has
* Prepares the @p SolutionTransfer for
* pure refinement. It
* stores the dof indices of each cell.
- * After calling this function
+ * After calling this function
* only calling the @p refine_interpolate
* functions is allowed.
*/
* coarsenend) grid.
*/
void prepare_for_coarsening_and_refinement (const std::vector<VECTOR> &all_in);
-
+
/**
* Same as previous function
* but for only one discrete function
* to be interpolated.
*/
void prepare_for_coarsening_and_refinement (const VECTOR &in);
-
+
/**
* This function
* interpolates the discrete function @p in,
*/
void refine_interpolate (const VECTOR &in,
VECTOR &out) const;
-
+
/**
* This function
* interpolates the discrete functions
*
* Calling this function is
* allowed only if first
- * Triangulation@p ::prepare_coarsening_and_refinement,
+ * Triangulation::prepare_coarsening_and_refinement,
* second
* @p SolutionTransfer::prepare_for_coarsening_and_refinement,
* an then third
- * Triangulation@p ::execute_coarsening_and_refinement
+ * Triangulation::execute_coarsening_and_refinement
* are called before. Multiple
* calling of this function is
* NOT allowed. Interpolating
*/
void interpolate (const std::vector<VECTOR>&all_in,
std::vector<VECTOR> &all_out) const;
-
+
/**
* Same as the previous function.
* It interpolates only one function.
* It assumes the vectors having the
* right sizes (i.e. <tt>in.size()==n_dofs_old</tt>,
* <tt>out.size()==n_dofs_refined</tt>)
- *
+ *
* Multiple calling of this function is
* NOT allowed. Interpolating
* several functions can be performed
* Exception
*/
DeclException0(ExcAlreadyPrepForCoarseAndRef);
-
+
/**
* Exception
*/
DeclException0(ExcTriaPrepCoarseningNotCalledBefore);
-
+
/**
* Exception
*/
* to work with.
*/
SmartPointer<const DH,SolutionTransfer<dim,VECTOR,DH> > dof_handler;
-
+
/**
* Stores the number of DoFs before the
* refinement and/or coarsening.
struct Pointerstruct {
Pointerstruct();
unsigned int memory_consumption () const;
-
+
std::vector<unsigned int> *indices_ptr;
std::vector<Vector<typename VECTOR::value_type> > *dof_values_ptr;
};
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2009, 2010 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* for i=0 to n_timesteps-1
* initialize timestep i for this sweep, e.g. for setting up
* data structures, creating temporary files, etc.
- *
+ *
* for i=0 to n_timesteps-1
* prepare timestep i for loop 0
* for i=0 to n_timesteps-1
* void TimeDependent_Wave<dim>::run_sweep (const unsigned int sweep_no)
* {
* start_sweep (sweep_no);
- *
+ *
* solve_primal_problem ();
- *
+ *
* if (compute_dual_problem)
* solve_dual_problem ();
- *
+ *
* postprocess ();
- *
+ *
* if (sweep_no != number_of_sweeps-1)
* refine_grids ();
- *
+ *
* write_statistics ();
- *
+ *
* end_sweep ();
* };
- *
- *
- *
+ *
+ *
+ *
* template <int dim>
- * void WaveProblem<dim>::run ()
+ * void WaveProblem<dim>::run ()
* {
* for (unsigned int sweep=0; sweep<number_of_sweeps; ++sweep)
* timestep_manager.run_sweep (sweep);
* };
* @endverbatim
* Here, @p timestep_manager is an object of type TimeDependent_Wave(), which
- * is a class derived from TimeDependent. @p start_sweep,
+ * is a class derived from TimeDependent. @p start_sweep,
* @p solve_primal_problem, @p solve_dual_problem, @p postprocess and @p end_sweep
- * are functions inherited from this class. They all do a loop over all
+ * are functions inherited from this class. They all do a loop over all
* timesteps within this object and call the respective function on each of
* these objects. For example, here are two of the functions as they are
* implemented by the library:
* @verbatim
- * void TimeDependent::start_sweep (const unsigned int s)
+ * void TimeDependent::start_sweep (const unsigned int s)
* {
* sweep_no = s;
- *
+ *
* // reset the number each
* // time step has, since some time
* // steps might have been added since
* timesteps[step]->set_timestep_no (step);
* timesteps[step]->set_sweep_no (sweep_no);
* };
- *
+ *
* for (unsigned int step=0; step<timesteps.size(); ++step)
* timesteps[step]->start_sweep ();
* };
- *
- *
+ *
+ *
* void
- * TimeDependent::solve_primal_problem ()
+ * TimeDependent::solve_primal_problem ()
* {
* do_loop (mem_fun(&TimeStepBase::init_for_primal_problem),
* mem_fun(&TimeStepBase::solve_primal_problem),
* invoked (@p solve_primal_problem, @p solve_dual_problem, @p postprocess,
* @p refine_grids and @p write_statistics all have this form, where the
* latter two give functions of the derived timestep class, rather than
- * from the base class). The function TimeStepBase@p ::init_for_primal_problem
+ * from the base class). The function TimeStepBase::init_for_primal_problem
* and the respective ones for the other operations defined by that class
* are only used to store the type of operation which the loop presently
* performed will do.
* for (unsigned int step=0; step<n_timesteps; ++step)
* init_function (static_cast<typename InitFunctionObject::argument_type>
* (timesteps[step]));
- *
+ *
* // wake up the first few time levels
* for (int step=-timestepping_data.look_ahead; step<0; ++step)
* for (int look_ahead=0; look_ahead<=timestepping_data.look_ahead; ++look_ahead)
* timesteps[step+look_ahead]->wake_up(look_ahead);
- *
- *
+ *
+ *
* for (unsigned int step=0; step<n_timesteps; ++step)
* {
* // first thing: wake up the
* for (unsigned int look_ahead=0;
* look_ahead<=timestepping_data.look_ahead; ++look_ahead)
* timesteps[step+look_ahead]->wake_up(look_ahead);
- *
- *
+ *
+ *
* // actually do the work
* loop_function (static_cast<typename LoopFunctionObject::argument_type>
* (timesteps[step]));
- *
+ *
* // let the timesteps behind sleep
* for (unsigned int look_back=0; look_back<=timestepping_data.look_back; ++look_back)
* timesteps[step-look_back]->sleep(look_back);
* };
- *
+ *
* // make the last few timesteps sleep
* for (int step=n_timesteps; n_timesteps+timestepping_data.look_back; ++step)
* for (int look_back=0; look_back<=timestepping_data.look_back; ++look_back)
*
* @author Wolfgang Bangerth, 1999
*/
-class TimeDependent
+class TimeDependent
{
public:
/**
*/
TimeSteppingData (const unsigned int look_ahead,
const unsigned int look_back);
-
+
/**
* This denotes the number of timesteps
* the timestepping algorithm needs to
* time step.
*/
const unsigned int look_ahead;
-
+
/**
* This is the opposite variable to the
* above one. It denotes the number of
* used in the previous sweep.
*/
void delete_timestep (const unsigned int position);
-
+
/**
* Solve the primal problem; uses the
* functions @p init_for_primal_problem
* object given to the constructor.
*/
void postprocess ();
-
+
/**
* Do a loop over all timesteps, call the
* @p init_function at the beginning and
*
* Instead of using the above form, you can
* equally well use
- * <tt>bind2nd(mem_fun1(&X::unary_function), arg)</tt>
+ * <tt>bind2nd(mem_fun1(&X::unary_function), arg)</tt>
* which lets the @p do_loop
* function call teh given function with
* the specified parameter. Note that you
* of this object.
*/
unsigned int memory_consumption () const;
-
+
/**
* Exception.
*/
DeclException0 (ExcInvalidPosition);
-
+
protected:
/**
* Vector holding pointers to the time
* do. See the documentation of this struct
* for more information.
*/
- const TimeSteppingData timestepping_data_postprocess;
+ const TimeSteppingData timestepping_data_postprocess;
private:
dual_problem = 0x1,
postprocess = 0x2
};
-
+
/**
* Constructor. Does nothing here apart
* from setting the time.
* and the like.
*/
virtual void end_sweep ();
-
+
/**
* Before the primal problem is
* solved on each time level, this
* are added or deleted.
*/
unsigned int get_timestep_no () const;
-
+
/**
* Compute the time difference to the
* last time step. If this timestep is
* above.
*/
double get_forward_timestep () const;
-
+
/**
* Determine an estimate for the
* memory consumption (in bytes)
* Exception
*/
DeclException0 (ExcCantComputeTimestep);
-
+
protected:
- /**
+ /**
* Pointer to the previous time step object
* in the list.
*/
const TimeStepBase *previous_timestep;
- /**
+ /**
* Pointer to the next time step object
* in the list.
*/
* functions.
*/
unsigned int next_action;
-
+
private:
/**
* Reset the pointer to the previous time
* really usable.
*/
Flags ();
-
+
/**
* Constructor; see the different
* fields for a description of the
Flags (const bool delete_and_rebuild_tria,
const unsigned int wakeup_level_to_build_grid,
const unsigned int sleep_level_to_delete_grid);
-
+
/**
* This flag determines whether
* the @p sleep and
* has no meaning.
*/
const unsigned int wakeup_level_to_build_grid;
-
+
/**
* This is the opposite flag to
* the one above: it determines
* the grid shall be deleted.
*/
const unsigned int sleep_level_to_delete_grid;
-
+
/**
* Exception
*/
* information.
*/
typedef std::vector<std::vector<std::pair<unsigned int, double> > > CorrectionRelaxations;
-
+
/**
* Default values for the relaxations:
* no relaxations.
*/
static CorrectionRelaxations default_correction_relaxations;
-
+
/**
* Constructor. The default
* values are chosen such that
const unsigned int cell_number_correction_steps = 0,
const bool mirror_flags_to_previous_grid = false,
const bool adapt_grids = false);
-
+
/**
* Maximum level of a cell in
* the triangulation of a time
* refinements.
*/
const unsigned int max_refinement_level;
-
+
/**
* First sweep to perform cell
* number correction steps on;
* this number of cells.
*/
const unsigned int min_cells_for_correction;
-
+
/**
* Fraction by which the number
* of cells on a time level may
* second: bottom deviation).
*/
const double cell_number_corridor_top;
-
+
/**
* @ref cell_number_corridor_top
*/
* correction step.
*/
const std::vector<std::vector<std::pair<unsigned int,double> > > correction_relaxations;
-
+
/**
* Number of iterations to be
* performed to adjust the
* previous one.
*/
const bool adapt_grids;
-
+
/**
* Exception
*/
* @p get_tria_refinement_criteria.
*/
template <int dim>
- struct RefinementData
+ struct RefinementData
{
/**
* Constructor
*/
RefinementData (const double refinement_threshold,
const double coarsening_threshold=0);
-
+
/**
* Threshold for refinement:
* cells having a larger value
* higher than this threshold).
*/
const double refinement_threshold;
-
+
/**
* Same threshold for
* coarsening: cells with a
typedef typename TimeStepBase_Tria_Flags::Flags<dim> Flags;
typedef typename TimeStepBase_Tria_Flags::RefinementFlags<dim> RefinementFlags;
typedef typename TimeStepBase_Tria_Flags::RefinementData<dim> RefinementData;
-
-
+
+
/**
* Extension of the enum in the base
* class denoting the next action to be
* terminal.
*/
TimeStepBase_Tria ();
-
+
/**
* Constructor. Takes a coarse
* grid from which the grids on this
* be overloaded.
*/
virtual void init_for_refinement ();
-
+
/**
* Virtual function that should fill
* the vector with the refinement
* the time step management object.
*/
virtual void get_tria_refinement_criteria (Vector<float> &criteria) const = 0;
-
+
/**
* The refinement
* flags of the triangulation are stored
DeclException0 (ExcGridNotDeleted);
protected:
-
+
/**
* Triangulation used at this
* time level. Since this is
* management object.
*/
SmartPointer<const Triangulation<dim, dim>,TimeStepBase_Tria<dim> > coarse_grid;
-
+
/**
* Some flags about how this time level
* shall behave. See the documentation
// some code at one place when it needs
// to be changed at a second place, here
// for the other direction, also.
-
+
const unsigned int n_timesteps = timesteps.size();
// initialize the time steps for
(&*timesteps[n_timesteps-step-1]));
break;
};
-
+
// let the timesteps behind sleep
for (unsigned int look_back=0;
look_back<=timestepping_data.look_back; ++look_back)
represented by the <code>Accessor</code> class. Usually, you will not use the
actual class names spelled out directly, but employ one of the typedefs
provided by the container classes, such as <code>typename
-Triangulation@<dim@>::cell_iterator</code>. Before going into this, let us
+Triangulation::cell_iterator</code>. Before going into this, let us
first discuss the concept of iterators, before delving into what the accessors
do.
<th>cell_iterator type</th>
<th>function call</th>
</tr>
-
+
<tr>
<th>Triangulation</th>
- <td>typename Triangulation@<dim,spacedim@>::cell_iterator</td>
+ <td>typename Triangulation::cell_iterator</td>
<td>triangulation.begin()</td>
</tr>
<tr>
<th>DoFHandler</th>
- <td>typename DoFHandler@<dim,spacedim@>::cell_iterator</td>
+ <td>typename DoFHandler::cell_iterator</td>
<td>dof_handler.begin()</td>
</tr>
<tr>
<th>hp::DoFHandler</th>
- <td>typename hp::DoFHandler@<dim,spacedim@>::cell_iterator</td>
+ <td>typename hp::DoFHandler::cell_iterator</td>
<td>hp_dof_handler.begin()</td>
</tr>
<tr>
<th>MGDoFHandler</th>
- <td>typename MGDoFHandler@<dim,spacedim@>::cell_iterator</td>
+ <td>typename MGDoFHandler::cell_iterator</td>
<td>mg_dof_handler.begin()</td>
</tr>
</table>
<th>face_iterator type</th>
<th>function call</th>
</tr>
-
+
<tr>
<th>Triangulation</th>
- <td>typename Triangulation@<dim,spacedim@>::face_iterator</td>
+ <td>typename Triangulation::face_iterator</td>
<td>triangulation.begin_face()</td>
</tr>
<tr>
<th>DoFHandler</th>
- <td>typename DoFHandler@<dim,spacedim@>::face_iterator</td>
+ <td>typename DoFHandler::face_iterator</td>
<td>dof_handler.begin_face()</td>
</tr>
<tr>
<th>hp::DoFHandler</th>
- <td>typename hp::DoFHandler@<dim,spacedim@>::face_iterator</td>
+ <td>typename hp::DoFHandler::face_iterator</td>
<td>hp_dof_handler.begin_face()</td>
</tr>
<tr>
<th>MGDoFHandler</th>
- <td>typename MGDoFHandler@<dim,spacedim@>::face_iterator</td>
+ <td>typename MGDoFHandler::face_iterator</td>
<td>mg_dof_handler.begin_face()</td>
</tr>
</table>
<th>cell_iterator type</th>
<th>function call</th>
</tr>
-
+
<tr>
<th>Triangulation</th>
- <td>typename Triangulation@<dim,spacedim@>::active_cell_iterator</td>
+ <td>typename Triangulation::active_cell_iterator</td>
<td>triangulation.begin_active()</td>
</tr>
<tr>
<th>DoFHandler</th>
- <td>typename DoFHandler@<dim,spacedim@>::active_cell_iterator</td>
+ <td>typename DoFHandler::active_cell_iterator</td>
<td>dof_handler.begin_active()</td>
</tr>
<tr>
<th>hp::DoFHandler</th>
- <td>typename hp::DoFHandler@<dim,spacedim@>::active_cell_iterator</td>
+ <td>typename hp::DoFHandler::active_cell_iterator</td>
<td>hp_dof_handler.begin_active()</td>
</tr>
<tr>
<th>MGDoFHandler</th>
- <td>typename MGDoFHandler@<dim,spacedim@>::active_cell_iterator</td>
+ <td>typename MGDoFHandler::active_cell_iterator</td>
<td>mg_dof_handler.begin_active()</td>
</tr>
</table>
<th>face_iterator type</th>
<th>function call</th>
</tr>
-
+
<tr>
<th>Triangulation</th>
- <td>typename Triangulation@<dim,spacedim@>::active_face_iterator</td>
+ <td>typename Triangulation::active_face_iterator</td>
<td>triangulation.begin_active_face()</td>
</tr>
<tr>
<th>DoFHandler</th>
- <td>typename DoFHandler@<dim,spacedim@>::active_face_iterator</td>
+ <td>typename DoFHandler::active_face_iterator</td>
<td>dof_handler.begin_active_face()</td>
</tr>
<tr>
<th>hp::DoFHandler</th>
- <td>typename hp::DoFHandler@<dim,spacedim@>::active_face_iterator</td>
+ <td>typename hp::DoFHandler::active_face_iterator</td>
<td>hp_dof_handler.begin_active_face()</td>
</tr>
<tr>
<th>MGDoFHandler</th>
- <td>typename MGDoFHandler@<dim,spacedim@>::active_face_iterator</td>
+ <td>typename MGDoFHandler::active_face_iterator</td>
<td>mg_dof_handler.begin_active_face()</td>
</tr>
</table>
internal::MGDoFHandler::Iterators<2,spacedim>, and
internal::MGDoFHandler::Iterators<3,spacedim> classes for MGDoFHandler
iterators.
-
-
+
+
@section IteratorAccessorInternals Iterator and accessor internals
Iterators, being like pointers, act as if they pointed to an actual
New: There is now a new RefinementCase
GeometryInfo::isotropic_refinement representing the isotropic
refinement case in <code>dim</code> dimensions,
- i.e. <tt>GeometryInfo<1>::isotropic_refinement=RefinementCase::cut_x</tt>,
- <tt>GeometryInfo<2>::isotropic_refinement=RefinementCase::cut_xy</tt>
+ i.e. GeometryInfo<1>::isotropic_refinement = RefinementCase::cut_x,
+ GeometryInfo<2>::isotropic_refinement = RefinementCase::cut_xy
and
- <tt>GeometryInfo<3>::isotropic_refinement=RefinementCase::cut_xyz</tt>.
+ GeometryInfo<3>::isotropic_refinement = RefinementCase::cut_xyz.
<br>
(RH 2008/07/08 as of branch_anisotropic at 2005/03/03)
</p>
/* $Id$ */
/* */
-/* Copyright (C) 2002, 2003, 2004, 2006, 2007, 2008, 2009 by the deal.II authors */
+/* Copyright (C) 2002, 2003, 2004, 2006, 2007, 2008, 2009, 2010 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
// for 1d. If the program is to be extended
// to higher space dimensions, so has to be
// this class.
-class InitializationValues : public Function<1>
+class InitializationValues : public Function<1>
{
public:
InitializationValues () : Function<1>() {}
-
+
virtual double value (const Point<1> &p,
const unsigned int component = 0) const;
};
// and return it, unless it is less than
// zero, in which case we take zero.
double InitializationValues::value (const Point<1> &p,
- const unsigned int) const
+ const unsigned int) const
{
const double base = std::pow(p(0), 1./3.);
const double random = 2.*rand()/RAND_MAX-1;
// step, as well as the present approximation
// of the solution.
template <int dim>
-class MinimizationProblem
+class MinimizationProblem
{
public:
MinimizationProblem (const unsigned int run_number);
void run ();
-
+
private:
void initialize_solution ();
void setup_system_on_mesh ();
static double energy (const DoFHandler<dim> &dof_handler,
const Vector<double> &function);
-
+
const unsigned int run_number;
-
+
Triangulation<dim> triangulation;
FE_Q<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix hanging_node_constraints;
-
+
SparsityPattern sparsity_pattern;
SparseMatrix<double> matrix;
// from erroneously compiling this function
// for other space dimensions, then.
template <>
-void MinimizationProblem<1>::initialize_solution ()
+void MinimizationProblem<1>::initialize_solution ()
{
// The first part is to assign the correct
// size to the vector, and use library
// the appropriate quadrature formula is
// the one we have chosen here.
QGauss<dim> quadrature_formula(4);
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
fe_values.shape_value(i,q_point)
*
(x_minus_u3 *
- (u_prime *
+ (u_prime *
fe_values.shape_grad(i,q_point))
-
(u_prime*u_prime) * u * u *
*
fe_values.JxW(q_point));
}
-
+
// After summing up all the
// contributions, we have to transfer
// them to the global objects. This is
matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
-
+
residual(local_dof_indices[i]) += cell_rhs(i);
}
}
// <code>alpha-dalpha</code> along the search
// direction:
const double dalpha = (alpha != 0 ? alpha/100 : 0.01);
-
+
tmp = present_solution;
tmp.add (alpha+dalpha, update);
const double f_a_plus = energy (dof_handler, tmp);
tmp = present_solution;
tmp.add (alpha+step_length, update);
const double e = energy (dof_handler, tmp);
-
+
if (e >= f_a)
step_length /= 2;
else
// the solution of a linear system before.
template <int dim>
void MinimizationProblem<dim>::do_step ()
-{
+{
assemble_step ();
Vector<double> update (present_solution.size());
SolverControl solver_control (residual.size(),
1e-2*residual.l2_norm());
SolverCG<> solver (solver_control);
-
+
PreconditionSSOR<> preconditioner;
preconditioner.initialize(matrix);
void MinimizationProblem<1>::refine_grid ()
{
const unsigned int dim = 1;
-
+
Vector<float> error_indicators (triangulation.n_active_cells());
// Then define the quadrature formula, and
5*(x-u*u*u)*u_doubleprime
+
2*u_prime*(1-3*u*u*u_prime)));
-
+
cell_residual_norm += (local_residual_value * local_residual_value *
fe_values.JxW(q));
}
DoFHandler<dim>::cell_iterator right_neighbor = cell->neighbor(1);
while (right_neighbor->has_children())
right_neighbor = right_neighbor->child(0);
-
+
neighbor_fe_values.reinit (right_neighbor);
neighbor_fe_values.get_function_grads (present_solution, local_gradients);
std::pow(u_prime_right,5));
error_indicators(cell_index) += right_jump * right_jump *
cell->diameter();
- }
- }
+ }
+ }
// Now we have all the refinement
// indicators computed, and want to refine
// situations, the library will silently
// also have to refine the neighbor cell
// once. It does so by calling the
- // <code>Triangulation@<dim@>::prepare_coarsening_and_refinement</code>
+ // <code>Triangulation::prepare_coarsening_and_refinement</code>
// function before actually doing the
// refinement and coarsening. This function
// flags a set of additional cells for
hanging_node_constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
hanging_node_constraints);
- hanging_node_constraints.close ();
+ hanging_node_constraints.close ();
hanging_node_constraints.distribute (present_solution);
// This is wasteful, since we create a
// <code>ConstraintMatrix</code> object that will be
// degree six, so a 4-point Gauss formula
// is appropriate:
QGauss<dim> quadrature_formula(4);
- FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
+ FEValues<dim> fe_values (dof_handler.get_fe(), quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
// converged somewhere, we output the
// results.
template <int dim>
-void MinimizationProblem<dim>::run ()
+void MinimizationProblem<dim>::run ()
{
GridGenerator::hyper_cube (triangulation, 0., 1.);
triangulation.refine_global (4);
initialize_solution ();
double last_energy = energy (dof_handler, present_solution);
-
+
while (true)
{
setup_system_on_mesh ();
-
+
for (unsigned int iteration=0; iteration<5; ++iteration)
do_step ();
refine_grid ();
}
-
+
output_results ();
-
+
std::cout << std::endl;
}
// time, all these runs are actually
// different, although it may seem that they
// are independent of each other.
-int main ()
+int main ()
{
try
{
for (unsigned int realization=0; realization<n_realizations; ++realization)
{
std::cout << "Realization " << realization << ":" << std::endl;
-
+
MinimizationProblem<1> minimization_problem_1d (realization);
minimization_problem_1d.run ();
}
<< std::endl;
return 1;
}
- catch (...)
+ catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
// linear system, we generate a
// preconditioner for the
// velocity-velocity matrix, i.e.,
- // <code>block(0,0)</code> in the system
- // matrix. As mentioned above, this
- // depends on the spatial
- // dimension. Since the two classes
- // described by the
- // <code>InnerPreconditioner@<dim@> ::
- // type</code> typedef have the same
- // interface, we do not have to do
- // anything different whether we want to
- // use a sparse direct solver or an ILU:
+ // <code>block(0,0)</code> in the
+ // system matrix. As mentioned
+ // above, this depends on the
+ // spatial dimension. Since the two
+ // classes described by the
+ // <code>InnerPreconditioner::type</code>
+ // typedef have the same interface,
+ // we do not have to do anything
+ // different whether we want to use
+ // a sparse direct solver or an
+ // ILU:
std::cout << " Computing preconditioner..." << std::endl << std::flush;
A_preconditioner
// introduction, the inverse is generated
// with the help of an inner preconditioner
// of type
- // <code>InnerPreconditioner@<dim@>::type</code>.
+ // <code>InnerPreconditioner::type</code>.
template <int dim>
void StokesProblem<dim>::solve ()
{
/* $Id$ */
/* Version: $Name: $ */
/* */
-/* Copyright (C) 2006, 2007, 2008, 2009 by the deal.II authors and Yaqi Wang */
+/* Copyright (C) 2006, 2007, 2008, 2009, 2010 by the deal.II authors and Yaqi Wang */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
{
public:
MaterialData (const unsigned int n_groups);
-
+
double get_diffusion_coefficient (const unsigned int group,
const unsigned int material_id) const;
double get_removal_XS (const unsigned int group,
const unsigned int material_id) const;
double get_fission_spectrum (const unsigned int group,
const unsigned int material_id) const;
-
+
private:
const unsigned int n_groups;
const unsigned int n_materials;
for (unsigned int group_2=0; group_2<n_groups; ++ group_2)
sigma_s[m][group_1][group_2] = 0.0;
}
-
-
+
+
diffusion[5][1] = 0.2;
-
+
sigma_r[4][0] = 0.026;
sigma_r[5][0] = 0.051;
sigma_r[6][0] = 0.026;
sigma_r[7][0] = 0.050;
-
+
sigma_r[0][1] = 0.100;
sigma_r[1][1] = 0.200;
sigma_r[2][1] = 0.250;
sigma_r[5][1] = 0.040;
sigma_r[6][1] = 0.020;
sigma_r[7][1] = 0.800;
-
+
nu_sigma_f[0][0] = 0.0050;
nu_sigma_f[1][0] = 0.0075;
nu_sigma_f[2][0] = 0.0075;
nu_sigma_f[5][0] = 0.000;
nu_sigma_f[6][0] = 1e-7;
nu_sigma_f[7][0] = 0.00;
-
+
nu_sigma_f[0][1] = 0.125;
nu_sigma_f[1][1] = 0.300;
nu_sigma_f[2][1] = 0.375;
nu_sigma_f[5][1] = 0.000;
nu_sigma_f[6][1] = 3e-6;
nu_sigma_f[7][1] = 0.00;
-
+
sigma_s[0][0][1] = 0.020;
sigma_s[1][0][1] = 0.015;
sigma_s[2][0][1] = 0.015;
default:
Assert (false,
ExcMessage ("Presently, only data for 2 groups is implemented"));
- }
+ }
}
MaterialData::get_diffusion_coefficient (const unsigned int group,
const unsigned int material_id) const
{
- Assert (group < n_groups,
+ Assert (group < n_groups,
ExcIndexRange (group, 0, n_groups));
- Assert (material_id < n_materials,
+ Assert (material_id < n_materials,
ExcIndexRange (material_id, 0, n_materials));
return diffusion[material_id][group];
MaterialData::get_removal_XS (const unsigned int group,
const unsigned int material_id) const
{
- Assert (group < n_groups,
+ Assert (group < n_groups,
ExcIndexRange (group, 0, n_groups));
- Assert (material_id < n_materials,
+ Assert (material_id < n_materials,
ExcIndexRange (material_id, 0, n_materials));
return sigma_r[material_id][group];
MaterialData::get_fission_XS (const unsigned int group,
const unsigned int material_id) const
{
- Assert (group < n_groups,
+ Assert (group < n_groups,
ExcIndexRange (group, 0, n_groups));
- Assert (material_id < n_materials,
+ Assert (material_id < n_materials,
ExcIndexRange (material_id, 0, n_materials));
return nu_sigma_f[material_id][group];
const unsigned int group_2,
const unsigned int material_id) const
{
- Assert (group_1 < n_groups,
+ Assert (group_1 < n_groups,
ExcIndexRange (group_1, 0, n_groups));
- Assert (group_2 < n_groups,
+ Assert (group_2 < n_groups,
ExcIndexRange (group_2, 0, n_groups));
- Assert (material_id < n_materials,
+ Assert (material_id < n_materials,
ExcIndexRange (material_id, 0, n_materials));
return sigma_s[material_id][group_1][group_2];
double
MaterialData::get_fission_spectrum (const unsigned int group,
- const unsigned int material_id) const
+ const unsigned int material_id) const
{
- Assert (group < n_groups,
+ Assert (group < n_groups,
ExcIndexRange (group, 0, n_groups));
- Assert (material_id < n_materials,
+ Assert (material_id < n_materials,
ExcIndexRange (material_id, 0, n_materials));
return chi[material_id][group];
double
MaterialData::get_fission_dist_XS (const unsigned int group_1,
const unsigned int group_2,
- const unsigned int material_id) const
+ const unsigned int material_id) const
{
return (get_fission_spectrum(group_1, material_id) *
get_fission_XS(group_2, material_id));
// include right hand side terms
// involving $\phi_g$. What one
// typically does in such situations
- // is to iterate: compute
+ // is to iterate: compute
// @f{eqnarray*}
// -\nabla \cdot(D_g(x) \nabla \phi^{(n)}_g(x))
// &+&
unsigned int n_active_cells () const;
unsigned int n_dofs () const;
-
+
// Then there are functions that
// assemble the linear system for
// each iteration and the present
// problem involving the neutron
// flux equations:
void solve ();
-
+
double get_fission_source () const;
-
+
void output_results (const unsigned int cycle) const;
-
+
void estimate_errors (Vector<float> &error_indicators) const;
-
+
void refine_grid (const Vector<float> &error_indicators,
const double refine_threshold,
const double coarsen_threshold);
// by the present guess of the
// eigenvalue we are looking for:
public:
-
+
Vector<double> solution;
Vector<double> solution_old;
// store them along with all the
// other data of this class:
private:
-
+
const unsigned int group;
const MaterialData &material_data;
-
+
Triangulation<dim> triangulation;
const FiniteElement<dim> &fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
-
+
Vector<double> system_rhs;
-
+
std::map<unsigned int,double> boundary_values;
ConstraintMatrix hanging_node_constraints;
-
+
// @sect5{Private member functionss}
//
EnergyGroup<dim>::setup_linear_system ()
{
const unsigned int n_dofs = dof_handler.n_dofs();
-
+
hanging_node_constraints.clear ();
DoFTools::make_hanging_node_constraints (dof_handler,
hanging_node_constraints);
hanging_node_constraints.close ();
system_matrix.clear ();
-
+
sparsity_pattern.reinit (n_dofs, n_dofs,
dof_handler.max_couplings_between_dofs());
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
// first clear the boundary value
// map):
boundary_values.clear();
-
+
for (unsigned int i=0; i<dim; ++i)
VectorTools::interpolate_boundary_values (dof_handler,
2*i+1,
// information is available.
template <int dim>
void
-EnergyGroup<dim>::assemble_system_matrix ()
+EnergyGroup<dim>::assemble_system_matrix ()
{
const QGauss<dim> quadrature_formula(fe.degree + 1);
-
- FEValues<dim> fe_values (fe, quadrature_formula,
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_JxW_values);
-
+
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
-
+
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
-
+
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
-
+
for (; cell!=endc; ++cell)
{
cell_matrix = 0;
-
+
fe_values.reinit (cell);
const double diffusion_coefficient
= material_data.get_diffusion_coefficient (group, cell->material_id());
const double removal_XS
= material_data.get_removal_XS (group,cell->material_id());
-
+
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
fe_values.shape_value(j,q_point))
*
fe_values.JxW(q_point));
-
+
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell; ++i)
local_dof_indices[j],
cell_matrix(i,j));
}
-
+
hanging_node_constraints.condense (system_matrix);
}
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
-
- FEValues<dim> fe_values (fe, quadrature_formula,
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_quadrature_points |
update_JxW_values);
-
+
Vector<double> cell_rhs (dofs_per_cell);
std::vector<double> extraneous_source_values (n_q_points);
std::vector<double> solution_old_values (n_q_points);
-
+
std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
-
+
for (; cell!=endc; ++cell)
{
cell_rhs = 0;
-
+
fe_values.reinit (cell);
-
+
const double fission_dist_XS
= material_data.get_fission_dist_XS (group, group, cell->material_id());
extraneous_source_values);
fe_values.get_function_values (solution_old, solution_old_values);
-
+
cell->get_dof_indices (local_dof_indices);
-
+
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
cell_rhs(i) += ((extraneous_source_values[q_point]
solution_old_values[q_point]) *
fe_values.shape_value(i,q_point) *
fe_values.JxW(q_point));
-
+
for (unsigned int i=0; i<dofs_per_cell; ++i)
- system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
{
if (group == g_prime.group)
return;
-
- const std::list<std::pair<typename DoFHandler<dim>::cell_iterator,
+
+ const std::list<std::pair<typename DoFHandler<dim>::cell_iterator,
typename DoFHandler<dim>::cell_iterator> >
cell_list
- = GridTools::get_finest_common_cells (dof_handler,
+ = GridTools::get_finest_common_cells (dof_handler,
g_prime.dof_handler);
-
- typename std::list<std::pair<typename DoFHandler<dim>::cell_iterator,
+
+ typename std::list<std::pair<typename DoFHandler<dim>::cell_iterator,
typename DoFHandler<dim>::cell_iterator> >
- ::const_iterator
+ ::const_iterator
cell_iter = cell_list.begin();
-
+
for (; cell_iter!=cell_list.end(); ++cell_iter)
{
FullMatrix<double> unit_matrix (fe.dofs_per_cell);
{
const QGauss<dim> quadrature_formula (fe.degree+1);
const unsigned int n_q_points = quadrature_formula.size();
-
- FEValues<dim> fe_values (fe, quadrature_formula,
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_JxW_values);
if (cell_g->level() > cell_g_prime->level())
fe_values.shape_value(i,q_point) *
fe_values.shape_value(j,q_point) *
fe_values.JxW(q_point));
- }
+ }
// Now we have all the
// interpolation (prolongation)
Vector<double> g_prime_old_values (fe.dofs_per_cell);
cell_g_prime->get_dof_values (g_prime.solution_old, g_prime_old_values);
cell_g_prime->get_dof_values (g_prime.solution, g_prime_new_values);
-
+
Vector<double> cell_rhs (fe.dofs_per_cell);
Vector<double> tmp (fe.dofs_per_cell);
local_mass_matrix_g.vmult (tmp, g_prime_new_values);
prolongation_matrix.Tvmult_add (cell_rhs, tmp);
- }
+ }
std::vector<unsigned int> local_dof_indices (fe.dofs_per_cell);
cell_g->get_dof_indices (local_dof_indices);
-
+
for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
{
const QGauss<dim> quadrature_formula (fe.degree + 1);
const unsigned int n_q_points = quadrature_formula.size();
-
- FEValues<dim> fe_values (fe, quadrature_formula,
+
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_JxW_values);
-
+
std::vector<double> solution_values (n_q_points);
double fission_source = 0;
-
+
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
= material_data.get_fission_XS(group, cell->material_id());
fe_values.get_function_values (solution, solution_values);
-
+
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
fission_source += (fission_XS *
solution_values[q_point] *
fe_values.JxW(q_point));
}
-
+
return fission_source;
}
// hand side vector.
template <int dim>
void
-EnergyGroup<dim>::solve ()
+EnergyGroup<dim>::solve ()
{
hanging_node_constraints.condense (system_rhs);
MatrixTools::apply_boundary_values (boundary_values,
system_matrix,
solution,
system_rhs);
-
+
SolverControl solver_control (system_matrix.m(),
1e-12*system_rhs.l2_norm());
SolverCG<> cg (solver_control);
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
-
+
cg.solve (system_matrix, solution, system_rhs, preconditioner);
hanging_node_constraints.distribute (solution);
void EnergyGroup<dim>::estimate_errors (Vector<float> &error_indicators) const
{
KellyErrorEstimator<dim>::estimate (dof_handler,
- QGauss<dim-1> (fe.degree + 1),
+ QGauss<dim-1> (fe.degree + 1),
typename FunctionMap<dim>::type(),
solution,
error_indicators);
typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
-
+
for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
if (error_indicators(cell_index) > refine_threshold)
cell->set_refine_flag ();
else if (error_indicators(cell_index) < coarsen_threshold)
cell->set_coarsen_flag ();
-
+
SolutionTransfer<dim> soltrans(dof_handler);
triangulation.prepare_coarsening_and_refinement();
Utilities::int_to_string(cycle,1) +
".eps";
std::ofstream output (filename.c_str());
-
+
GridOut grid_out;
grid_out.write_eps (triangulation, output);
}
"." +
Utilities::int_to_string(cycle,1) +
".gmv";
-
+
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
-
+
std::ofstream output (filename.c_str());
data_out.write_gmv (output);
}
// parameters from a ParameterHandler
// object.
template <int dim>
-class NeutronDiffusionProblem
+class NeutronDiffusionProblem
{
public:
class Parameters
{
public:
Parameters ();
-
+
static void declare_parameters (ParameterHandler &prm);
void get_parameters (ParameterHandler &prm);
unsigned int n_refinement_cycles;
unsigned int fe_degree;
-
+
double convergence_tolerance;
};
-
-
+
+
NeutronDiffusionProblem (const Parameters ¶meters);
~NeutronDiffusionProblem ();
void run ();
-
+
private:
// @sect5{Private member functions}
-
+
// There are not that many member
// functions in this class since
// most of the functionality has
-
+
// @sect4{Implementation of the <code>NeutronDiffusionProblem</code> class}
// Now for the
template <int dim>
-NeutronDiffusionProblem<dim>::~NeutronDiffusionProblem ()
+NeutronDiffusionProblem<dim>::~NeutronDiffusionProblem ()
{
for (unsigned int group=0; group<energy_groups.size(); ++group)
delete energy_groups[group];
:
Point<dim> (assemblies_x*rods_per_assembly_x*pin_pitch_x,
assemblies_y*rods_per_assembly_y*pin_pitch_y,
- assemblies_z*assembly_height));
+ assemblies_z*assembly_height));
std::vector<unsigned int> n_subdivisions;
n_subdivisions.push_back (assemblies_x*rods_per_assembly_x);
n_subdivisions.push_back (assemblies_y*rods_per_assembly_y);
if (dim >= 3)
n_subdivisions.push_back (assemblies_z);
-
+
Triangulation<dim> coarse_grid;
GridGenerator::subdivided_hyper_rectangle (coarse_grid,
n_subdivisions,
upper_right,
true);
-
+
// The second part of the function
// deals with material numbers of
// pin cells of each type of
{ 1, 1, 1, 1, 1, 5, 1, 1, 5, 1, 1, 5, 1, 1, 1, 1, 1 },
{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 }
- },
+ },
{
{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 2, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 3, 3, 3, 2 },
{ 2, 3, 3, 5, 3, 4, 4, 4, 4, 4, 4, 4, 3, 5, 3, 3, 2 },
{ 2, 3, 3, 3, 3, 5, 3, 3, 5, 3, 3, 5, 3, 3, 3, 3, 2 },
- { 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2 },
+ { 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 2 },
{ 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2 }
- },
+ },
{
{ 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
{ 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 },
// the 'UX' and 'PX' assemblies:
const unsigned int core[assemblies_x][assemblies_y][assemblies_z]
= {{{0}, {2}}, {{2}, {0}}};
-
+
// We are now in a position to
// actually set material IDs for
// each cell. To this end, we loop
++cell)
{
const Point<dim> cell_center = cell->center();
-
+
const unsigned int tmp_x = int(cell_center[0]/pin_pitch_x);
const unsigned int ax = tmp_x/rods_per_assembly_x;
const unsigned int cx = tmp_x - ax * rods_per_assembly_x;
Assert (cx < rods_per_assembly_x, ExcInternalError());
Assert (cy < rods_per_assembly_y, ExcInternalError());
-
+
cell->set_material_id(assembly_materials[core[ax][ay][az]][cx][cy] - 1);
}
// spawned returns a double. We can then later
// ask this object for the returned value
// (when doing so, the
- // Threads::Thread@<double@>::return_value
+ // Threads::Thread::return_value
// function first waits for the thread to
// finish if it hasn't done so already).
//
for (unsigned int group=0; group<parameters.n_groups; ++group)
threads.push_back (Threads::new_thread (&EnergyGroup<dim>::get_fission_source,
*energy_groups[group]));
-
+
double fission_source = 0;
for (unsigned int group=0; group<parameters.n_groups; ++group)
fission_source += threads[group].return_value ();
BlockVector<float> group_error_indicators(n_cells);
- {
+ {
Threads::ThreadGroup<> threads;
for (unsigned int group=0; group<parameters.n_groups; ++group)
threads += Threads::new_thread (&EnergyGroup<dim>::estimate_errors,
group_error_indicators.block(group));
threads.join_all ();
}
-
+
const float max_error = group_error_indicators.linfty_norm();
const float refine_threshold = 0.3*max_error;
const float coarsen_threshold = 0.01*max_error;
- {
+ {
Threads::ThreadGroup<> threads;
for (unsigned int group=0; group<parameters.n_groups; ++group)
threads += Threads::new_thread (&EnergyGroup<dim>::refine_grid,
- *energy_groups[group],
+ *energy_groups[group],
group_error_indicators.block(group),
refine_threshold,
coarsen_threshold);
// there is actually not much to
// comment on:
template <int dim>
-void NeutronDiffusionProblem<dim>::run ()
+void NeutronDiffusionProblem<dim>::run ()
{
std::cout << std::setprecision (12) << std::fixed;
Timer timer;
timer.start ();
-
+
for (unsigned int cycle=0; cycle<parameters.n_refinement_cycles; ++cycle)
{
std::cout << "Cycle " << cycle << ':' << std::endl;
-
+
if (cycle == 0)
initialize_problem();
else
refine_grid ();
for (unsigned int group=0; group<parameters.n_groups; ++group)
energy_groups[group]->solution *= k_eff;
- }
+ }
for (unsigned int group=0; group<parameters.n_groups; ++group)
energy_groups[group]->setup_linear_system ();
std::cout << energy_groups[group]->n_dofs()
<< ' ';
std::cout << std::endl << std::endl;
-
-
+
+
Threads::ThreadGroup<> threads;
for (unsigned int group=0; group<parameters.n_groups; ++group)
threads += Threads::new_thread
for (unsigned int group=0; group<parameters.n_groups; ++group)
{
energy_groups[group]->assemble_ingroup_rhs (ZeroFunction<dim>());
-
+
for (unsigned int bgroup=0; bgroup<parameters.n_groups; ++bgroup)
energy_groups[group]->assemble_cross_group_rhs (*energy_groups[bgroup]);
-
+
energy_groups[group]->solve ();
}
while((error > parameters.convergence_tolerance)
&&
(iteration < 500));
-
+
for (unsigned int group=0; group<parameters.n_groups; ++group)
energy_groups[group]->output_results (cycle);
<< ", time=" << timer()
<< std::endl;
-
+
std::cout << std::endl << std::endl;
}
}
// of type
// <code>NeutronDiffusionProblem</code>
// for computation of the eigenvalue:
-int main (int argc, char ** argv)
+int main (int argc, char ** argv)
{
const unsigned int dim = 2;
-
+
try
{
deallog.depth_console (0);
else
filename = argv[1];
-
+
ParameterHandler parameter_handler;
-
+
NeutronDiffusionProblem<dim>::Parameters parameters;
parameters.declare_parameters (parameter_handler);
parameter_handler.read_input (filename);
-
+
parameters.get_parameters (parameter_handler);
-
+
NeutronDiffusionProblem<dim> neutron_diffusion_problem (parameters);
neutron_diffusion_problem.run ();
}
return 1;
}
- catch (...)
+ catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
An error occurred in line <273> of file <step-5.cc> in function
void Coefficient<dim>::value_list(const std::vector<Point<dim>, std::allocator<Point<dim> > >&, std::vector<double, std::allocator<double> >&, unsigned int)
const [with int dim = 2]
-The violated condition was:
+The violated condition was:
values.size() == points.size()
The name and call sequence of the exception was:
ExcDimensionMismatch (values.size(), points.size())
-Additional Information:
+Additional Information:
Dimension 1 not equal to 2
Stacktrace:
An error occurred in line <273> of file <step-5.cc> in function
void Coefficient<dim>::value_list(const std::vector<Point<dim>, std::allocator<Point<dim> > >&, std::vector<double, std::allocator<double> >&, unsigned int)
const [with int dim = 2]
-The violated condition was:
+The violated condition was:
values.size() == points.size()
The name and call sequence of the exception was:
ExcDimensionMismatch (values.size(), points.size())
-Additional Information:
+Additional Information:
Dimension 1 not equal to 2
Stacktrace:
So now we know pretty exactly where the error happened and what went
wrong. What we don't know yet is how exactly we got there. The
stacktrace at the bottom actually tells us what happened: the problem
-happened in
+happened in
<code>Coefficient::value_list</code> (stackframe 0) and that it was
called from <code>main()</code> (stackframe 1). In realistic programs,
there would be many more functions in between these two. For example,
we might have made the mistake in the <code>assemble_system</code>
function, in which case stack frame 1 would be
-<code>LaplaceProblem<2>::assemble_system</code>, stack frame 2
-would be <code>LaplaceProblem<2>::run</code>, and stack frame 3
+<code>LaplaceProblem::assemble_system</code>, stack frame 2
+would be <code>LaplaceProblem::run</code>, and stack frame 3
would be <code>main()</code> — you get the idea.
// clear when we discuss its
// implementation.
template <int dim>
-class LaplaceProblem
+class LaplaceProblem
{
public:
LaplaceProblem ();
~LaplaceProblem ();
void run ();
-
+
private:
void setup_system ();
void assemble_system ();
// from step-5:
template <int dim>
-class Coefficient : public Function<dim>
+class Coefficient : public Function<dim>
{
public:
Coefficient () : Function<dim>() {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
-
+
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component = 0) const;
template <int dim>
double Coefficient<dim>::value (const Point<dim> &p,
- const unsigned int) const
+ const unsigned int) const
{
if (p.square() < 0.5*0.5)
return 20;
template <int dim>
void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
- const unsigned int component) const
+ const unsigned int component) const
{
const unsigned int n_points = points.size();
- Assert (values.size() == n_points,
+ Assert (values.size() == n_points,
ExcDimensionMismatch (values.size(), n_points));
-
- Assert (component == 0,
+
+ Assert (component == 0,
ExcIndexRange (component, 0, 1));
-
+
for (unsigned int i=0; i<n_points; ++i)
{
if (points[i].square() < 0.5*0.5)
// destructor, to the end of the
// results section of this example.
template <int dim>
-LaplaceProblem<dim>::~LaplaceProblem ()
+LaplaceProblem<dim>::~LaplaceProblem ()
{
dof_handler.clear ();
}
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
-
+
// After setting up all the degrees
// of freedoms, here are now the
// differences compared to step-5,
// sparsity pattern immediately.
CompressedSparsityPattern c_sparsity(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern (dof_handler, c_sparsity);
-
+
// The constrained hanging nodes
// will later be eliminated from
// the linear system of
// code and nothing that you have to
// worry about.
template <int dim>
-void LaplaceProblem<dim>::assemble_system ()
-{
+void LaplaceProblem<dim>::assemble_system ()
+{
const QGauss<dim> quadrature_formula(3);
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
coefficient.value_list (fe_values.get_quadrature_points(),
coefficient_values);
-
+
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
-
+
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
// find at the end of this function:
template <int dim>
-void LaplaceProblem<dim>::solve ()
+void LaplaceProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
// corresponding Neumann values. This
// information is represented by an
// object of type
- // <code>FunctionMap@<dim@>::type</code> that is
+ // <code>FunctionMap::type</code> that is
// essentially a map from boundary
// indicators to function objects
// describing Neumann boundary values
std::string filename = "grid-";
filename += ('0' + cycle);
filename += ".eps";
-
+
std::ofstream output (filename.c_str());
GridOut grid_out;
// The rest of the loop looks as
// before:
template <int dim>
-void LaplaceProblem<dim>::run ()
+void LaplaceProblem<dim>::run ()
{
for (unsigned int cycle=0; cycle<8; ++cycle)
{
}
else
refine_grid ();
-
+
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
-
+
assemble_system ();
solve ();
output_results (cycle);
// factor of four.
DataOutBase::EpsFlags eps_flags;
eps_flags.z_scaling = 4;
-
+
DataOut<dim> data_out;
data_out.set_flags (eps_flags);
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
-
+
std::ofstream output ("final-solution.eps");
data_out.write_eps (output);
}
// actually encodes the functionality
// particular to the present
// application.
-int main ()
+int main ()
{
// The general idea behind the
// use the <code>what</code> member function
// to get a string which describes
// the reason why the exception was
- // thrown.
+ // thrown.
//
// The deal.II exception classes
// are all derived from the
// can't do anything at all. We
// then simply print an error
// message and exit.
- catch (...)
+ catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
/* $Id$ */
/* */
-/* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008 by the deal.II authors */
+/* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
// below in the constructor of this
// class.
template <int dim>
-class ElasticProblem
+class ElasticProblem
{
public:
ElasticProblem ();
~ElasticProblem ();
void run ();
-
+
private:
void setup_system ();
void assemble_system ();
// the scalar function used in
// previous programs.
template <int dim>
-class RightHandSide : public Function<dim>
+class RightHandSide : public Function<dim>
{
public:
RightHandSide ();
-
+
// The next change is that we
// want a replacement for the
// <code>value</code> function of the
template <int dim>
inline
void RightHandSide<dim>::vector_value (const Point<dim> &p,
- Vector<double> &values) const
+ Vector<double> &values) const
{
- Assert (values.size() == dim,
+ Assert (values.size() == dim,
ExcDimensionMismatch (values.size(), dim));
Assert (dim >= 2, ExcNotImplemented());
-
+
// The rest of the function
// implements computing force
// values. We will use a constant
Point<dim> point_1, point_2;
point_1(0) = 0.5;
point_2(0) = -0.5;
-
+
// If now the point <code>p</code> is in a
// circle (sphere) of radius 0.2
// around one of these points, then
values(0) = 1;
else
values(0) = 0;
-
+
// Likewise, if <code>p</code> is in the
// vicinity of the origin, then set
// the y-force to 1, otherwise to
if (p.square() < 0.2*0.2)
values(1) = 1;
else
- values(1) = 0;
+ values(1) = 0;
}
// make some things simpler below.
template <int dim>
void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
- std::vector<Vector<double> > &value_list) const
+ std::vector<Vector<double> > &value_list) const
{
- Assert (value_list.size() == points.size(),
+ Assert (value_list.size() == points.size(),
ExcDimensionMismatch (value_list.size(), points.size()));
const unsigned int n_points = points.size();
//
// We can prevent this situation by
// calling
- // <code>RightHandSide@<dim@>::vector_valued</code>
+ // <code>RightHandSide::vector_value</code>
// on each point in the input
// list. Note that by giving the
// full name of the function,
// The destructor, on the other hand,
// is exactly as in step-6:
template <int dim>
-ElasticProblem<dim>::~ElasticProblem ()
+ElasticProblem<dim>::~ElasticProblem ()
{
dof_handler.clear ();
}
// this is not explicit knowledge we
// need to care about:
template <int dim>
-void ElasticProblem<dim>::assemble_system ()
-{
+void ElasticProblem<dim>::assemble_system ()
+{
QGauss<dim> quadrature_formula(2);
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
cell_rhs = 0;
fe_values.reinit (cell);
-
+
// Next we get the values of
// the coefficients at the
// quadrature points. Likewise
right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
rhs_values);
-
+
// Then assemble the entries of
// the local stiffness matrix
// and right hand side
// contributions:
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
- const unsigned int
+ const unsigned int
component_i = fe.system_to_component_index(i).first;
-
- for (unsigned int j=0; j<dofs_per_cell; ++j)
+
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
{
- const unsigned int
+ const unsigned int
component_j = fe.system_to_component_index(j).first;
-
+
for (unsigned int q_point=0; q_point<n_q_points;
++q_point)
{
- cell_matrix(i,j)
- +=
+ cell_matrix(i,j)
+ +=
// The first term
// is (lambda d_i
// u_i, d_j v_j)
// introduction:
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
- const unsigned int
+ const unsigned int
component_i = fe.system_to_component_index(i).first;
-
+
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
cell_rhs(i) += fe_values.shape_value(i,q_point) *
rhs_values[q_point](component_i) *
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
-
+
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
}
// the system indeed is. Therefore,
// we need not change anything.
template <int dim>
-void ElasticProblem<dim>::solve ()
+void ElasticProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverCG<> cg (solver_control);
std::string filename = "solution-";
filename += ('0' + cycle);
Assert (cycle < 10, ExcInternalError());
-
+
filename += ".gmv";
std::ofstream output (filename.c_str());
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
-
+
// As said above, we need a
// different name for each
solution_names.push_back ("displacement");
break;
case 2:
- solution_names.push_back ("x_displacement");
+ solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
break;
case 3:
- solution_names.push_back ("x_displacement");
+ solution_names.push_back ("x_displacement");
solution_names.push_back ("y_displacement");
solution_names.push_back ("z_displacement");
break;
default:
Assert (false, ExcNotImplemented());
}
-
+
// After setting up the names for
// the different components of the
// solution vector, we can add the
// cells without actually refining
// the grid.)
template <int dim>
-void ElasticProblem<dim>::run ()
+void ElasticProblem<dim>::run ()
{
for (unsigned int cycle=0; cycle<8; ++cycle)
{
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
-
+
assemble_system ();
solve ();
output_results (cycle);
// The main function is again exactly
// like in step-6 (apart from the
// changed class names, of course).
-int main ()
+int main ()
{
try
{
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
-
+
return 1;
}
- catch (...)
+ catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
/* $Id$ */
/* */
-/* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008 by the deal.II authors */
+/* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007, 2008, 2010 by the deal.II authors */
/* */
/* This file is subject to QPL and may not be distributed */
/* without copyright and license information. Please refer */
// previous examples, so we again
// only comment on the differences.
template <int dim>
-class AdvectionProblem
+class AdvectionProblem
{
public:
AdvectionProblem ();
~AdvectionProblem ();
void run ();
-
+
private:
void setup_system ();
// The next function will be used
void assemble_system ();
void assemble_system_interval (const typename DoFHandler<dim>::active_cell_iterator &begin,
const typename DoFHandler<dim>::active_cell_iterator &end);
-
+
// The following functions again
// are as in previous examples,
// as are the subsequent
{
public:
AdvectionField () : TensorFunction<1,dim> () {}
-
+
virtual Tensor<1,dim> value (const Point<dim> &p) const;
-
+
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<Tensor<1,dim> > &values) const;
// assertion as in this case can
// eliminate many of these problems.
template <int dim>
-Tensor<1,dim>
-AdvectionField<dim>::value (const Point<dim> &p) const
+Tensor<1,dim>
+AdvectionField<dim>::value (const Point<dim> &p) const
{
Point<dim> value;
value[0] = 2;
template <int dim>
void
AdvectionField<dim>::value_list (const std::vector<Point<dim> > &points,
- std::vector<Tensor<1,dim> > &values) const
+ std::vector<Tensor<1,dim> > &values) const
{
Assert (values.size() == points.size(),
ExcDimensionMismatch (values.size(), points.size()));
-
+
for (unsigned int i=0; i<points.size(); ++i)
values[i] = AdvectionField<dim>::value (points[i]);
}
{
public:
RightHandSide () : Function<dim>() {}
-
+
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
-
+
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component = 0) const;
-
+
private:
static const Point<dim> center_point;
};
template <int dim>
double
RightHandSide<dim>::value (const Point<dim> &p,
- const unsigned int component) const
+ const unsigned int component) const
{
Assert (component == 0, ExcIndexRange (component, 0, 1));
const double diameter = 0.1;
void
RightHandSide<dim>::value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
- const unsigned int component) const
+ const unsigned int component) const
{
Assert (values.size() == points.size(),
ExcDimensionMismatch (values.size(), points.size()));
-
+
for (unsigned int i=0; i<points.size(); ++i)
values[i] = RightHandSide<dim>::value (points[i], component);
}
virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
-
+
virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component = 0) const;
template <int dim>
double
BoundaryValues<dim>::value (const Point<dim> &p,
- const unsigned int component) const
+ const unsigned int component) const
{
Assert (component == 0, ExcIndexRange (component, 0, 1));
void
BoundaryValues<dim>::value_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
- const unsigned int component) const
+ const unsigned int component) const
{
Assert (values.size() == points.size(),
ExcDimensionMismatch (values.size(), points.size()));
-
+
for (unsigned int i=0; i<points.size(); ++i)
values[i] = BoundaryValues<dim>::value (points[i], component);
}
static void estimate_interval (const DoFHandler<dim> &dof,
const Vector<double> &solution,
const IndexInterval &index_interval,
- Vector<float> &error_per_cell);
+ Vector<float> &error_per_cell);
};
template <int dim>
-AdvectionProblem<dim>::~AdvectionProblem ()
+AdvectionProblem<dim>::~AdvectionProblem ()
{
dof_handler.clear ();
}
// same size and assembling on each
// of these chunks in parallel.
template <int dim>
-void AdvectionProblem<dim>::assemble_system ()
+void AdvectionProblem<dim>::assemble_system ()
{
// First, we want to find out how
// many threads shall assemble the
// only for a range of cell
// iterators, but for iterators in
// general, so you could use it for
- // <code>std::vector@<T@>::iterator</code> or
+ // <code>std::vector::iterator</code> or
// usual pointers as well.
//
// The function returns a vector of
// an alias.
typedef typename DoFHandler<dim>::active_cell_iterator active_cell_iterator;
std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
- thread_ranges
+ thread_ranges
= Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
dof_handler.end (),
n_threads);
// multi-threading, then no threads
// can run in parallel and the
// function returns immediately.
- threads.join_all ();
+ threads.join_all ();
// After the matrix has been
}
-
+
// Now, this is the function that
// does the actual work. It is not
// very different from the
void
AdvectionProblem<dim>::
assemble_system_interval (const typename DoFHandler<dim>::active_cell_iterator &begin,
- const typename DoFHandler<dim>::active_cell_iterator &end)
+ const typename DoFHandler<dim>::active_cell_iterator &end)
{
// First of all, we will need some
// objects that describe boundary
const AdvectionField<dim> advection_field;
const RightHandSide<dim> right_hand_side;
const BoundaryValues<dim> boundary_values;
-
+
// Next we need quadrature formula
// for the cell terms, but also for
// the integral over the inflow
// direction are sufficient.
QGauss<dim> quadrature_formula(2);
QGauss<dim-1> face_quadrature_formula(2);
-
+
// Finally, we need objects of type
// <code>FEValues</code> and
// <code>FEFaceValues</code>. For the cell
// integrals, we don't need the
// gradients, but rather the normal
// vectors to the cells.
- FEValues<dim> fe_values (fe, quadrature_formula,
+ FEValues<dim> fe_values (fe, quadrature_formula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
// FEFaceValues object
// for the present face:
fe_face_values.reinit (cell, face);
-
+
// For the quadrature
// points at hand, we ask
// for the values of the
face_boundary_values);
advection_field.value_list (fe_face_values.get_quadrature_points(),
face_advection_directions);
-
+
// Now loop over all
// quadrature points and
// see whether it is on
fe_face_values.shape_value(i,q_point) *
fe_face_values.shape_value(j,q_point) *
fe_face_values.JxW(q_point));
-
+
cell_rhs(i) -= (face_advection_directions[q_point] *
fe_face_values.normal_vector(q_point) *
face_boundary_values[q_point] *
fe_face_values.JxW(q_point));
};
};
-
+
// Now go on by transferring
// the local contributions to
system_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i,j));
-
+
system_rhs(local_dof_indices[i]) += cell_rhs(i);
};
assembler_lock.release ();
// method. As preconditioner, we use
// the Jacobi method.
template <int dim>
-void AdvectionProblem<dim>::solve ()
+void AdvectionProblem<dim>::solve ()
{
SolverControl solver_control (1000, 1e-12);
SolverBicgstab<> bicgstab (solver_control);
std::string filename = "grid-";
filename += ('0' + cycle);
Assert (cycle < 10, ExcInternalError());
-
+
filename += ".eps";
std::ofstream output (filename.c_str());
// ... as is the main loop (setup --
// solve -- refine)
template <int dim>
-void AdvectionProblem<dim>::run ()
+void AdvectionProblem<dim>::run ()
{
for (unsigned int cycle=0; cycle<6; ++cycle)
{
{
refine_grid ();
};
-
+
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs()
<< std::endl;
-
+
assemble_system ();
solve ();
output_results (cycle);
data_out.attach_dof_handler (dof_handler);
data_out.add_data_vector (solution, "solution");
data_out.build_patches ();
-
+
std::ofstream output ("final-solution.gmv");
data_out.write_gmv (output);
}
// except for delegating work to the
// other function:
template <int dim>
-void
+void
GradientEstimation::estimate (const DoFHandler<dim> &dof_handler,
const Vector<double> &solution,
Vector<float> &error_per_cell)
//
// Now for the details:
template <int dim>
-void
+void
GradientEstimation::estimate_interval (const DoFHandler<dim> &dof_handler,
const Vector<double> &solution,
const IndexInterval &index_interval,
FEValues<dim> fe_midpoint_value (dof_handler.get_fe(),
midpoint_rule,
update_values | update_quadrature_points);
-
+
// Then we need space foe the
// tensor <code>Y</code>, which is the sum
// of outer products of the
// element for this iterator.
Vector<float>::iterator
error_on_this_cell = error_per_cell.begin() + index_interval.first;
-
+
// Then we allocate a vector to
// hold iterators to all active
// abbreviation for the
// iterator to the face
// and the neighbor
- const typename DoFHandler<dim>::face_iterator
+ const typename DoFHandler<dim>::face_iterator
face = cell->face(face_no);
- const typename DoFHandler<dim>::cell_iterator
+ const typename DoFHandler<dim>::cell_iterator
neighbor = cell->neighbor(face_no);
// Then check whether the
neighbor_child = neighbor;
while (neighbor_child->has_children())
neighbor_child = neighbor_child->child (face_no==0 ? 1 : 0);
-
+
// As this used
// some
// non-trivial
// here.
Assert (neighbor_child->neighbor(face_no==0 ? 1 : 0)==cell,
ExcInternalError());
-
+
// If the check
// succeeded, we
// push the
std::vector<double> this_midpoint_value(1);
fe_midpoint_value.get_function_values (solution, this_midpoint_value);
-
+
// Now loop over all active neighbors
// and collect the data we
// neighbor cell:
const typename DoFHandler<dim>::active_cell_iterator
neighbor = *neighbor_ptr;
-
+
// Then get the center of
// the neighbor cell and
// the value of the finite
Point<dim> y = neighbor_center - this_center;
const double distance = std::sqrt(y.square());
y /= distance;
-
+
// Then add up the
// contribution of this
// cell to the Y matrix...
for (unsigned int i=0; i<dim; ++i)
for (unsigned int j=0; j<dim; ++j)
Y[i][j] += y[i] * y[j];
-
+
// ... and update the sum
// of difference quotients:
projected_gradient += (neighbor_midpoint_value[0] -
// right powers of the mesh
// width:
const Tensor<2,dim> Y_inverse = invert(Y);
-
+
Point<dim> gradient;
contract (gradient, Y_inverse, projected_gradient);
-
+
*error_on_this_cell = (std::pow(cell->diameter(),
1+1.0*dim/2) *
std::sqrt(gradient.square()));
// the only difference in the name of
// the main class that actually does
// the computation.
-int main ()
+int main ()
{
try
{
<< std::endl;
return 1;
}
- catch (...)
+ catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
* filter. Constrained dofs are
* transformed to local index space of
* the filter, and elements not present
- * in the IndexSet are ignored.
+ * in the IndexSet are ignored.
*
* This function provides an easy way to
* create a ConstraintMatrix for certain
* set_inhomogeneity().
*/
void add_lines (const IndexSet &lines);
-
+
/**
* Add an entry to a given
* line. The list of lines is
* Does the same as the function above
* but can treat
* non quadratic matrices.
- */
+ */
template <typename MatrixType>
void
distribute_local_to_global (const FullMatrix<double> &local_matrix,
* vector instead of a map and the
* consequences thereof, the same
* applies as what is said for
- * ConstraintMatrix@p ::lines.
+ * ConstraintMatrix::lines.
*/
Entries entries;
* Krylov space method inside the preconditioner.
*
* FGMRES needs two vectors in each iteration steps yielding a total
- * of <tt>2*AdditionalData::max_basis_size+1</tt> auxiliary vectors.
+ * of <tt>2 * SolverFGMRESAdditionalData::max_basis_size+1</tt>
+ * auxiliary vectors.
*
* Caveat: documentation of this class is not up to date. There are
* also a few parameters of GMRES we would like to introduce here.
const double s = vv.l2_norm();
h(inner_iteration+1) = s;
//TODO: s=0 is a lucky breakdown. Handle this somehow decently
-
+
vv *= 1./s;
/* Transformation into