* Read or write the data of this object for serialization
*/
template <class Archive>
- void serialize(Archive &ar, const unsigned int version) {
+ void serialize(Archive &ar, const unsigned int version)
+ {
ar &valid
- &h5_sol_filename
- &h5_mesh_filename
- &entry_time
- &num_nodes
- &num_cells
- &dimension
- &attribute_dims;
+ &h5_sol_filename
+ &h5_mesh_filename
+ &entry_time
+ &num_nodes
+ &num_cells
+ &dimension
+ &attribute_dims;
}
/// Get the XDMF content associated with this entry.
/**
* Override the standard function that returns the description of the error.
*/
- virtual const char* what() const throw();
+ virtual const char *what() const throw();
/**
* Get exception name.
// Fill the fields of the exception object
e.set_fields (file, line, function, cond, exc_name);
- switch(handling)
+ switch (handling)
{
case abort_on_exception:
dealii::deal_II_exceptions::internals::abort(e);
*/
#ifdef DEBUG
#define Assert(cond, exc) \
-{ \
- if (!(cond)) \
- ::dealii::deal_II_exceptions::internals:: \
- issue_error(::dealii::deal_II_exceptions::internals::abort_on_exception,\
- __FILE__, __LINE__, __PRETTY_FUNCTION__, #cond, #exc, exc); \
-}
+ { \
+ if (!(cond)) \
+ ::dealii::deal_II_exceptions::internals:: \
+ issue_error(::dealii::deal_II_exceptions::internals::abort_on_exception,\
+ __FILE__, __LINE__, __PRETTY_FUNCTION__, #cond, #exc, exc); \
+ }
#else
#define Assert(cond, exc) \
-{}
+ {}
#endif
*/
#ifdef DEBUG
#define AssertNothrow(cond, exc) \
-{ \
- if (!(cond)) \
- ::dealii::deal_II_exceptions::internals:: \
- issue_error( \
- ::dealii::deal_II_exceptions::internals::abort_nothrow_on_exception, \
- __FILE__, __LINE__, __PRETTY_FUNCTION__, #cond, #exc, exc); \
-}
+ { \
+ if (!(cond)) \
+ ::dealii::deal_II_exceptions::internals:: \
+ issue_error( \
+ ::dealii::deal_II_exceptions::internals::abort_nothrow_on_exception, \
+ __FILE__, __LINE__, __PRETTY_FUNCTION__, #cond, #exc, exc); \
+ }
#else
#define AssertNothrow(cond, exc) \
-{}
+ {}
#endif
*/
#ifdef HAVE_BUILTIN_EXPECT
#define AssertThrow(cond, exc) \
-{ \
- if (__builtin_expect(!(cond), false)) \
- ::dealii::deal_II_exceptions::internals:: \
- issue_error(::dealii::deal_II_exceptions::internals::throw_on_exception,\
- __FILE__, __LINE__, __PRETTY_FUNCTION__, #cond, #exc, exc); \
-}
+ { \
+ if (__builtin_expect(!(cond), false)) \
+ ::dealii::deal_II_exceptions::internals:: \
+ issue_error(::dealii::deal_II_exceptions::internals::throw_on_exception,\
+ __FILE__, __LINE__, __PRETTY_FUNCTION__, #cond, #exc, exc); \
+ }
#else /*ifdef HAVE_BUILTIN_EXPECT*/
#define AssertThrow(cond, exc) \
-{ \
- if (!(cond)) \
- ::dealii::deal_II_exceptions::internals:: \
- issue_error(::dealii::deal_II_exceptions::internals::throw_on_exception,\
- __FILE__, __LINE__, __PRETTY_FUNCTION__, #cond, #exc, exc); \
-}
+ { \
+ if (!(cond)) \
+ ::dealii::deal_II_exceptions::internals:: \
+ issue_error(::dealii::deal_II_exceptions::internals::throw_on_exception,\
+ __FILE__, __LINE__, __PRETTY_FUNCTION__, #cond, #exc, exc); \
+ }
#endif /*ifdef HAVE_BUILTIN_EXPECT*/
static void
alternating_form_at_vertices
#ifndef DEAL_II_CONSTEXPR_BUG
- (const Point<spacedim> (&vertices)[vertices_per_cell],
- Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell]);
+ (const Point<spacedim> (&vertices)[vertices_per_cell],
+ Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell]);
#else
- (const Point<spacedim> *vertices,
- Tensor<spacedim-dim,spacedim> *forms);
+ (const Point<spacedim> *vertices,
+ Tensor<spacedim-dim,spacedim> *forms);
#endif
/**
Assert (i < dim, ExcIndexRange(i, 0, dim));
static const RefinementCase options[dim] = { RefinementPossibilities<2>::cut_x,
- RefinementPossibilities<2>::cut_y };
+ RefinementPossibilities<2>::cut_y
+ };
return options[i];
}
Assert (i < dim, ExcIndexRange(i, 0, dim));
static const RefinementCase options[dim] = { RefinementPossibilities<3>::cut_x,
- RefinementPossibilities<3>::cut_y,
- RefinementPossibilities<3>::cut_z };
+ RefinementPossibilities<3>::cut_y,
+ RefinementPossibilities<3>::cut_z
+ };
return options[i];
}
return (p[0] <= 0.5 ?
(p[1] <= 0.5 ? 0 : 2) :
- (p[1] <= 0.5 ? 1 : 3));
+ (p[1] <= 0.5 ? 1 : 3));
}
return (p[0] <= 0.5 ?
(p[1] <= 0.5 ?
(p[2] <= 0.5 ? 0 : 4) :
- (p[2] <= 0.5 ? 2 : 6)) :
- (p[1] <= 0.5 ?
- (p[2] <= 0.5 ? 1 : 5) :
- (p[2] <= 0.5 ? 3 : 7)));
+ (p[2] <= 0.5 ? 2 : 6)) :
+ (p[1] <= 0.5 ?
+ (p[2] <= 0.5 ? 1 : 5) :
+ (p[2] <= 0.5 ? 3 : 7)));
}
// intersection of the conflict indices.
for (unsigned int i=0; i<partition_size; ++i)
for (unsigned int j=i+1; j<partition_size; ++j)
- // If the two iterators share indices then we increase the degree of the
- // vertices and create an ''edge'' in the graph.
- if (have_nonempty_intersection (conflict_indices[i], conflict_indices[j]))
- {
- ++degrees[i];
- ++degrees[j];
- graph[i].push_back(j);
- graph[j].push_back(i);
- }
+ // If the two iterators share indices then we increase the degree of the
+ // vertices and create an ''edge'' in the graph.
+ if (have_nonempty_intersection (conflict_indices[i], conflict_indices[j]))
+ {
+ ++degrees[i];
+ ++degrees[j];
+ graph[i].push_back(j);
+ graph[j].push_back(i);
+ }
// Sort the vertices by decreasing degree.
std::vector<int>::iterator degrees_it;
*/
Threads::ThreadLocalStorage<std_cxx1x::shared_ptr<std::ostringstream> > outstreams;
- template <typename T> friend LogStream &operator << (LogStream & log, const T &t);
+ template <typename T> friend LogStream &operator << (LogStream &log, const T &t);
};
*/
template <typename T>
inline
-LogStream & operator<< (LogStream &log, const T &t)
+LogStream &operator<< (LogStream &log, const T &t)
{
// print to the internal stringstream
log.get_stream() << t;
* three elements.
*/
template <typename I1, typename I2,
- typename I3, typename I4>
+ typename I3, typename I4>
static
void
apply (const F &f,
* module.
*/
template <typename InputIterator,
- typename OutputIterator,
- typename Predicate>
+ typename OutputIterator,
+ typename Predicate>
void transform (const InputIterator &begin_in,
const InputIterator &end_in,
OutputIterator out,
* module.
*/
template <typename InputIterator1,
- typename InputIterator2,
- typename OutputIterator,
- typename Predicate>
+ typename InputIterator2,
+ typename OutputIterator,
+ typename Predicate>
void transform (const InputIterator1 &begin_in1,
const InputIterator1 &end_in1,
InputIterator2 in2,
* module.
*/
template <typename InputIterator1,
- typename InputIterator2,
- typename InputIterator3,
- typename OutputIterator,
- typename Predicate>
+ typename InputIterator2,
+ typename InputIterator3,
+ typename OutputIterator,
+ typename Predicate>
void transform (const InputIterator1 &begin_in1,
const InputIterator1 &end_in1,
InputIterator2 in2,
* for each subrange.
*/
template <typename ResultType,
- typename Function>
+ typename Function>
struct ReductionOnSubranges
{
/**
* Advance a tuple of iterators by $n$.
*/
template <typename I1, typename I2,
- typename I3, typename I4>
+ typename I3, typename I4>
inline
void advance (std_cxx1x::tuple<I1,I2,I3, I4> &t,
const unsigned int n)
* Advance a tuple of iterators by 1.
*/
template <typename I1, typename I2,
- typename I3, typename I4>
+ typename I3, typename I4>
inline
void advance_by_one (std_cxx1x::tuple<I1,I2,I3,I4> &t)
{
* for each number of arguments.
*/
template <typename RT, typename ArgList,
- int length = std_cxx1x::tuple_size<ArgList>::value>
+
int length = std_cxx1x::tuple_size<ArgList>::value>
struct fun_ptr_helper;
* returned by the function called on
* the thread.
*/
- std_cxx1x::shared_ptr<return_value<RT> > ret_val;
+ std_cxx1x::shared_ptr<return_value<RT> > ret_val;
/**
* Start the thread and
* Copy constructor.
*/
Thread (const Thread<RT> &t)
- :
- thread_descriptor (t.thread_descriptor)
+ :
+ thread_descriptor (t.thread_descriptor)
{}
/**
* @deprecated Use new_thread() instead.
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3>,3>
spawn (RT (*fun_ptr)(Arg1,Arg2,Arg3)) DEAL_II_DEPRECATED;
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3>,3>
spawn (RT (*fun_ptr)(Arg1,Arg2,Arg3))
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3>,3>
spawn (C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3)) DEAL_II_DEPRECATED;
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3>,3>
spawn (C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3))
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3>,3>
spawn (const C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3) const) DEAL_II_DEPRECATED;
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3>,3>
spawn (const C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3) const)
* @deprecated Use new_thread() instead.
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4>,4>
spawn (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4)) DEAL_II_DEPRECATED;
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4>,4>
spawn (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4))
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4>,4>
spawn (C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4)) DEAL_II_DEPRECATED;
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4>,4>
spawn (C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4))
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4>,4>
spawn (const C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4) const) DEAL_II_DEPRECATED;
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4>,4>
spawn (const C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4) const)
* @deprecated Use new_thread() instead.
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5>,5>
spawn (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5)) DEAL_II_DEPRECATED;
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5>,5>
spawn (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5))
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5>,5>
spawn (C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5)) DEAL_II_DEPRECATED;
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5>,5>
spawn (C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5))
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5>,5>
spawn (const C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5) const) DEAL_II_DEPRECATED;
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5>,5>
spawn (const C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5) const)
* @deprecated Use new_thread() instead.
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5, Arg6>,6>
spawn (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6)) DEAL_II_DEPRECATED;
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5, Arg6>,6>
spawn (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6))
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5, Arg6>,6>
spawn (C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6)) DEAL_II_DEPRECATED;
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5, Arg6>,6>
spawn (C &c, RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6))
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
internal::fun_encapsulator<RT,
std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5, Arg6>,6>
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
internal::fun_encapsulator<RT,
std_cxx1x::tuple<Arg1, Arg2, Arg3, Arg4, Arg5, Arg6>,6>
* @deprecated Use new_thread() instead.
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6, Arg7>,7>
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6, Arg7>,7>
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6, Arg7>,7>
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6, Arg7>,7>
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
internal::fun_encapsulator<RT,
std_cxx1x::tuple<Arg1, Arg2, Arg3,
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
internal::fun_encapsulator<RT,
std_cxx1x::tuple<Arg1, Arg2, Arg3,
* @deprecated Use new_thread() instead.
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6,
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6,
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6,
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6,
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
internal::fun_encapsulator<RT,
std_cxx1x::tuple<Arg1, Arg2, Arg3,
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
internal::fun_encapsulator<RT,
std_cxx1x::tuple<Arg1, Arg2, Arg3,
* @deprecated Use new_thread() instead.
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6,
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6,
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6,
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
internal::fun_encapsulator<RT,std_cxx1x::tuple<Arg1, Arg2, Arg3,
Arg4, Arg5, Arg6,
* @deprecated Use new_thread() instead.
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
internal::fun_encapsulator<RT,
std_cxx1x::tuple<Arg1, Arg2, Arg3,
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
internal::fun_encapsulator<RT,
std_cxx1x::tuple<Arg1, Arg2, Arg3,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
Thread<RT>
new_thread (RT (*fun_ptr)(Arg1,Arg2,Arg3),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
Thread<RT>
new_thread (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
Thread<RT>
new_thread (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
Thread<RT>
new_thread (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
Thread<RT>
new_thread (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6,Arg7),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6,Arg7),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6,Arg7) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
Thread<RT>
new_thread (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
Thread<RT>
new_thread (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
Thread<RT>
new_thread (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* Copy constructor.
*/
Task (const Task<RT> &t)
- :
- task_descriptor (t.task_descriptor)
+ :
+ task_descriptor (t.task_descriptor)
{}
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
Task<RT>
new_task (RT (*fun_ptr)(Arg1,Arg2,Arg3),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3>
+ typename Arg1, typename Arg2, typename Arg3>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
Task<RT>
new_task (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3, typename Arg4>
+ typename Arg1, typename Arg2, typename Arg3, typename Arg4>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
Task<RT>
new_task (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
Task<RT>
new_task (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
Task<RT>
new_task (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6,Arg7),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6,Arg7),
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,Arg6,Arg7) const,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
Task<RT>
new_task (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
Task<RT>
new_task (RT (*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* @ingroup threads
*/
template <typename RT, typename C,
- typename Arg1, typename Arg2, typename Arg3,
- typename Arg4, typename Arg5, typename Arg6,
- typename Arg7, typename Arg8, typename Arg9>
+ typename Arg1, typename Arg2, typename Arg3,
+ typename Arg4, typename Arg5, typename Arg6,
+ typename Arg7, typename Arg8, typename Arg9>
inline
Task<RT>
new_task (RT (C::*fun_ptr)(Arg1,Arg2,Arg3,Arg4,Arg5,
* A class that creates a sequence of items from a range of iterators.
*/
template <typename Iterator,
- typename ScratchData,
- typename CopyData>
+ typename ScratchData,
+ typename CopyData>
class IteratorRangeToItemStream : public tbb::filter
{
public:
* run in parallel.
*/
template <typename Iterator,
- typename ScratchData,
- typename CopyData>
+ typename ScratchData,
+ typename CopyData>
class Worker : public tbb::filter
{
public:
* they are created.
*/
template <typename Iterator,
- typename ScratchData,
- typename CopyData>
+ typename ScratchData,
+ typename CopyData>
class Copier : public tbb::filter
{
public:
* A class that creates a sequence of items from a range of iterators.
*/
template <typename Iterator,
- typename ScratchData,
- typename CopyData>
+ typename ScratchData,
+ typename CopyData>
class IteratorRangeToItemStream : public tbb::filter
{
public:
* run in parallel.
*/
template <typename Iterator,
- typename ScratchData,
- typename CopyData>
+ typename ScratchData,
+ typename CopyData>
class WorkerAndCopier : public tbb::filter
{
public:
* are generated.
*/
template <typename Worker,
- typename Copier,
- typename Iterator,
- typename ScratchData,
- typename CopyData>
+ typename Copier,
+ typename Iterator,
+ typename ScratchData,
+ typename CopyData>
void
run (const Iterator &begin,
const typename identity<Iterator>::type &end,
* are generated.
*/
template <typename Worker,
- typename Copier,
- typename Iterator,
- typename ScratchData,
- typename CopyData>
+ typename Copier,
+ typename Iterator,
+ typename ScratchData,
+ typename CopyData>
void
run (const std::vector<std::vector<Iterator> > &colored_iterators,
Worker worker,
* are generated.
*/
template <typename MainClass,
- typename Iterator,
- typename ScratchData,
- typename CopyData>
+ typename Iterator,
+ typename ScratchData,
+ typename CopyData>
void
run (const Iterator &begin,
const typename identity<Iterator>::type &end,
}
//forward declaration of the data type for periodic face pairs
-namespace GridTools {template <typename CellIterator> struct PeriodicFacePair;}
+namespace GridTools
+{
+ template <typename CellIterator> struct PeriodicFacePair;
+}
namespace parallel
{
register_data_attach (const std::size_t size,
const std_cxx1x::function<void (const cell_iterator &,
const CellStatus,
- void *)> & pack_callback);
+ void *)> &pack_callback);
/**
* The given function is called for
notify_ready_to_unpack (const unsigned int offset,
const std_cxx1x::function<void (const cell_iterator &,
const CellStatus,
- const void *)> & unpack_callback);
+ const void *)> &unpack_callback);
/**
* Returns a permutation vector for the order the coarse
*/
void
add_periodicity
- (const std::vector<GridTools::PeriodicFacePair<cell_iterator> >&);
+ (const std::vector<GridTools::PeriodicFacePair<cell_iterator> > &);
*/
void
fill_vertices_with_ghost_neighbors
- (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- &vertices_with_ghost_neighbors);
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors);
template <int, int> friend class dealii::internal::DoFHandler::Policy::ParallelDistributed;
};
*/
void
fill_vertices_with_ghost_neighbors
- (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- &vertices_with_ghost_neighbors);
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors);
};
}
{
return dof_handler.faces->lines.
get_dof_index (dof_handler,
- obj_index,
+ obj_index,
fe_index,
local_index,
obj_level);
const dealii::internal::int2type<1> &)
{
return dof_handler.levels[obj_level]->fe_index_is_active(obj_index,
- fe_index);
+ fe_index);
}
const unsigned int obj_index,
const dealii::internal::int2type<1> &)
{
- // on a cell, the number of active elements is one
+ // on a cell, the number of active elements is one
return 1;
}
const unsigned int n,
const dealii::internal::int2type<1> &)
{
- Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
+ Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
return dof_handler.levels[obj_level]->active_fe_index (obj_index);
}
const dealii::internal::int2type<2> &)
{
return dof_handler.levels[obj_level]->fe_index_is_active(obj_index,
- fe_index);
+ fe_index);
}
const unsigned int obj_index,
const dealii::internal::int2type<2> &)
{
- // on a cell, the number of active elements is one
+ // on a cell, the number of active elements is one
return 1;
}
const unsigned int n,
const dealii::internal::int2type<2> &)
{
- Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
+ Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
return dof_handler.levels[obj_level]->active_fe_index (obj_index);
}
const dealii::internal::int2type<3> &)
{
return dof_handler.levels[obj_level]->fe_index_is_active(obj_index,
- fe_index);
+ fe_index);
}
const unsigned int obj_index,
const dealii::internal::int2type<3> &)
{
- // on a cell, the number of active elements is one
+ // on a cell, the number of active elements is one
return 1;
}
const unsigned int n,
const dealii::internal::int2type<3> &)
{
- Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
+ Assert (n==0, ExcMessage("On cells, there can only be one active FE index"));
return dof_handler.levels[obj_level]->active_fe_index (obj_index);
}
Assert (static_cast<unsigned int>(accessor.present_index)
<
accessor.dof_handler->levels[accessor.present_level]
- ->cell_cache_offsets.size(),
+ ->cell_cache_offsets.size(),
ExcInternalError());
Assert (accessor.dof_handler->levels[accessor.present_level]
->cell_cache_offsets[accessor.present_index]
std::vector<types::global_dof_index> dof_indices (dofs_per_cell);
static_cast<const dealii::DoFAccessor<dim,dealii::hp::DoFHandler<dim,spacedim>,lda> &>
- (accessor).get_dof_indices (dof_indices, accessor.active_fe_index());
+ (accessor).get_dof_indices (dof_indices, accessor.active_fe_index());
types::global_dof_index *next_dof_index
- = &accessor.dof_handler->levels[accessor.present_level]
- ->cell_dof_indices_cache[accessor.dof_handler->levels[accessor.present_level]
- ->cell_cache_offsets[accessor.present_index]];
+ = &accessor.dof_handler->levels[accessor.present_level]
+ ->cell_dof_indices_cache[accessor.dof_handler->levels[accessor.present_level]
+ ->cell_cache_offsets[accessor.present_index]];
for (unsigned int i=0; i<dofs_per_cell; ++i, ++next_dof_index)
*next_dof_index = dof_indices[i];
}
ExcMessage ("DoFHandler not initialized"));
return accessor.dof_handler->levels[accessor.level()]
- ->active_fe_index(accessor.present_index);
+ ->active_fe_index(accessor.present_index);
}
ExcMessage ("DoFHandler not initialized"));
accessor.dof_handler->levels[accessor.level()]
- ->set_active_fe_index (accessor.present_index, i);
+ ->set_active_fe_index (accessor.present_index, i);
}
template <int dim, int spacedim, bool lda, typename number,
- class OutputMatrix, typename OutputVector>
+ class OutputMatrix, typename OutputVector>
static
void
distribute_local_to_global (const DoFCellAccessor<dealii::DoFHandler<dim,spacedim>, lda> &accessor,
template <int dim, int spacedim, bool lda, typename number,
- class OutputMatrix, typename OutputVector>
+ class OutputMatrix, typename OutputVector>
static
void
distribute_local_to_global (const DoFCellAccessor<dealii::hp::DoFHandler<dim,spacedim>, lda> &accessor,
AssertDimension (dof_indices.size(), this->get_fe().dofs_per_cell);
const types::global_dof_index *cache
- = this->dof_handler->levels[this->present_level]
- ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
+ = this->dof_handler->levels[this->present_level]
+ ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
for (unsigned int i=0; i<this->get_fe().dofs_per_cell; ++i, ++cache)
dof_indices[i] = *cache;
}
typename DoFCellAccessor::ExcVectorDoesNotMatch());
const types::global_dof_index *cache
- = this->dof_handler->levels[this->present_level]
- ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
+ = this->dof_handler->levels[this->present_level]
+ ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
values.extract_subvector_to (cache,
cache + this->get_fe().dofs_per_cell,
== this->get_fe().dofs_per_cell,
typename DoFCellAccessor::ExcVectorDoesNotMatch());
Assert (values.size() == this->get_dof_handler().n_dofs(),
- typename DoFCellAccessor::ExcVectorDoesNotMatch());
+ typename DoFCellAccessor::ExcVectorDoesNotMatch());
const types::global_dof_index *cache
- = this->dof_handler->levels[this->present_level]
- ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
+ = this->dof_handler->levels[this->present_level]
+ ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
constraints.get_dof_values(values, *cache, local_values_begin,
local_values_end);
== this->get_fe().dofs_per_cell,
typename DoFCellAccessor::ExcVectorDoesNotMatch());
Assert (values.size() == this->get_dof_handler().n_dofs(),
- typename DoFCellAccessor::ExcVectorDoesNotMatch());
+ typename DoFCellAccessor::ExcVectorDoesNotMatch());
const types::global_dof_index *cache
- = this->dof_handler->levels[this->present_level]
- ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
+ = this->dof_handler->levels[this->present_level]
+ ->get_cell_cache_start (this->present_index, this->get_fe().dofs_per_cell);
for (unsigned int i=0; i<this->get_fe().dofs_per_cell; ++i, ++cache)
values(*cache) = local_values(i);
*
* Defined in dof_handler.cc.
*/
- template<int dim, int spacedim>
- std::string policy_to_string(const dealii::internal::DoFHandler::Policy::PolicyBase<dim,spacedim> & policy);
+ template<int dim, int spacedim>
+ std::string policy_to_string(const dealii::internal::DoFHandler::Policy::PolicyBase<dim,spacedim> &policy);
}
"the one that was associated with the DoFHandler previously stored."));
AssertThrow (policy_name == internal::policy_to_string(*policy),
ExcMessage (std::string ("The policy currently associated with this DoFHandler (")
- + internal::policy_to_string(*policy)
- +std::string(") does not match the one that was associated with the "
- "DoFHandler previously stored (")
- + policy_name
- + ")."));
+ + internal::policy_to_string(*policy)
+ +std::string(") does not match the one that was associated with the "
+ "DoFHandler previously stored (")
+ + policy_name
+ + ")."));
}
*/
const types::global_dof_index *
get_cell_cache_start (const unsigned int obj_index,
- const unsigned int dofs_per_cell) const;
+ const unsigned int dofs_per_cell) const;
/**
* Determine an estimate for the
inline
const types::global_dof_index *
DoFLevel<dim>::get_cell_cache_start (const unsigned int obj_index,
- const unsigned int dofs_per_cell) const
+ const unsigned int dofs_per_cell) const
{
Assert (obj_index*dofs_per_cell+dofs_per_cell
<=
template <class GridClass> class InterGridMap;
template <int dim, int spacedim> class Mapping;
-namespace GridTools {template <typename CellIterator> struct PeriodicFacePair;}
+namespace GridTools
+{
+ template <typename CellIterator> struct PeriodicFacePair;
+}
//TODO: map_support_points_to_dofs should generate a multimap, rather than just a map, since several dofs may be located at the same support point
template<typename DH>
void
make_periodicity_constraints
- (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
- &periodic_faces,
+ (const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
+ &periodic_faces,
dealii::ConstraintMatrix &constraint_matrix,
const ComponentMask &component_mask = ComponentMask());
* @p distribute function of your hanging node constraints object.
*/
template <int dim, int spacedim,
- template <int,int> class DH1,
- template <int,int> class DH2,
- class InVector, class OutVector>
+ template <int,int> class DH1,
+ template <int,int> class DH2,
+ class InVector, class OutVector>
void
interpolate (const DH1<dim,spacedim> &dof1,
const InVector &u1,
* DoF values at the discontinuities.
*/
template <int dim, int spacedim,
- template <int, int> class DH1,
- template <int, int> class DH2,
- class InVector, class OutVector>
+ template <int, int> class DH1,
+ template <int, int> class DH2,
+ class InVector, class OutVector>
void interpolate (const DH1<dim,spacedim> &dof1,
const InVector &u1,
const DH2<dim,spacedim> &dof2,
* type @p hp::DoFHandler.
*/
template <int dim,
- template <int> class DH,
- class InVector, class OutVector, int spacedim>
+ template <int> class DH,
+ class InVector, class OutVector, int spacedim>
void back_interpolate (const DH<dim> &dof1,
const InVector &u1,
const FiniteElement<dim,spacedim> &fe2,
*/
DeclException1 (ExcAccessToUninitializedField,
char *,
- << ("You are requesting information from an FEValues/FEFaceValues/FESubfaceValues "
- "object for which this kind of information has not been computed. What "
- "information these objects compute is determined by the update_* flags you "
- "pass to the constructor. Here, the operation you are attempting requires "
- "the <")
- << arg1
- << "> flag to be set, but it was apparently not specified upon construction.");
+ << ("You are requesting information from an FEValues/FEFaceValues/FESubfaceValues "
+ "object for which this kind of information has not been computed. What "
+ "information these objects compute is determined by the update_* flags you "
+ "pass to the constructor. Here, the operation you are attempting requires "
+ "the <")
+ << arg1
+ << "> flag to be set, but it was apparently not specified upon construction.");
/**
* @todo Document this
*
FEValuesBase<dim,spacedim>::get_quadrature_points () const
{
Assert (this->update_flags & update_quadrature_points,
- ExcAccessToUninitializedField("update_quadrature_points"));
+ ExcAccessToUninitializedField("update_quadrature_points"));
return this->quadrature_points;
}
FEValuesBase<dim,spacedim>::get_JxW_values () const
{
Assert (this->update_flags & update_JxW_values,
- ExcAccessToUninitializedField("update_JxW_values"));
+ ExcAccessToUninitializedField("update_JxW_values"));
return this->JxW_values;
}
FEValuesBase<dim,spacedim>::get_jacobians () const
{
Assert (this->update_flags & update_jacobians,
- ExcAccessToUninitializedField("update_jacobians"));
+ ExcAccessToUninitializedField("update_jacobians"));
return this->jacobians;
}
FEValuesBase<dim,spacedim>::get_jacobian_grads () const
{
Assert (this->update_flags & update_jacobian_grads,
- ExcAccessToUninitializedField("update_jacobians_grads"));
+ ExcAccessToUninitializedField("update_jacobians_grads"));
return this->jacobian_grads;
}
FEValuesBase<dim,spacedim>::get_inverse_jacobians () const
{
Assert (this->update_flags & update_inverse_jacobians,
- ExcAccessToUninitializedField("update_inverse_jacobians"));
+ ExcAccessToUninitializedField("update_inverse_jacobians"));
return this->inverse_jacobians;
}
FEValuesBase<dim,spacedim>::quadrature_point (const unsigned int i) const
{
Assert (this->update_flags & update_quadrature_points,
- ExcAccessToUninitializedField("update_quadrature_points"));
+ ExcAccessToUninitializedField("update_quadrature_points"));
Assert (i<this->quadrature_points.size(), ExcIndexRange(i, 0, this->quadrature_points.size()));
return this->quadrature_points[i];
FEValuesBase<dim,spacedim>::JxW (const unsigned int i) const
{
Assert (this->update_flags & update_JxW_values,
- ExcAccessToUninitializedField("update_JxW_values"));
+ ExcAccessToUninitializedField("update_JxW_values"));
Assert (i<this->JxW_values.size(), ExcIndexRange(i, 0, this->JxW_values.size()));
return this->JxW_values[i];
FEValuesBase<dim,spacedim>::jacobian (const unsigned int i) const
{
Assert (this->update_flags & update_jacobians,
- ExcAccessToUninitializedField("update_jacobians"));
+ ExcAccessToUninitializedField("update_jacobians"));
Assert (i<this->jacobians.size(), ExcIndexRange(i, 0, this->jacobians.size()));
return this->jacobians[i];
FEValuesBase<dim,spacedim>::jacobian_grad (const unsigned int i) const
{
Assert (this->update_flags & update_jacobian_grads,
- ExcAccessToUninitializedField("update_jacobians_grads"));
+ ExcAccessToUninitializedField("update_jacobians_grads"));
Assert (i<this->jacobian_grads.size(), ExcIndexRange(i, 0, this->jacobian_grads.size()));
return this->jacobian_grads[i];
FEValuesBase<dim,spacedim>::inverse_jacobian (const unsigned int i) const
{
Assert (this->update_flags & update_inverse_jacobians,
- ExcAccessToUninitializedField("update_inverse_jacobians"));
+ ExcAccessToUninitializedField("update_inverse_jacobians"));
Assert (i<this->inverse_jacobians.size(), ExcIndexRange(i, 0, this->inverse_jacobians.size()));
return this->inverse_jacobians[i];
*/
template<int dim>
void hyper_cube_with_cylindrical_hole (
- Triangulation<dim> &triangulation,
- const double inner_radius = .25,
- const double outer_radius = .5,
- const double L = .5,
- const unsigned int repetition = 1,
- const bool colorize = false);
+ Triangulation<dim> &triangulation,
+ const double inner_radius = .25,
+ const double outer_radius = .5,
+ const double L = .5,
+ const unsigned int repetition = 1,
+ const bool colorize = false);
/**
* Produce a ring of cells in 3D that is cut open, twisted and glued
void serialize(Archive &ar,
const unsigned int version);
};
-
+
/**
* Store boundary and material data. For
* example, in one dimension, this field
return *(pointer + 1 + local_index);
else
pointer += static_cast<types::global_dof_index>(
- dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<structdim>() + 1);
+ dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>() + 1);
}
}
return true;
else
pointer += static_cast<types::global_dof_index>(
- dof_handler.get_fe()[*pointer]
- .template n_dofs_per_object<structdim>()+1);
+ dof_handler.get_fe()[*pointer]
+ .template n_dofs_per_object<structdim>()+1);
}
}
* <tt>levels[]</tt> tree of
* the Triangulation objects.
*/
- std::vector<dealii::internal::hp::DoFLevel*> levels;
+ std::vector<dealii::internal::hp::DoFLevel *> levels;
/**
* Space to store the DoF
*/
void
set_active_fe_index (const unsigned int obj_index,
- const unsigned int fe_index);
+ const unsigned int fe_index);
/**
* Return a pointer to the beginning of the DoF indices cache
*/
const types::global_dof_index *
get_cell_cache_start (const unsigned int obj_index,
- const unsigned int dofs_per_cell) const;
+ const unsigned int dofs_per_cell) const;
/**
* Determine an estimate for the
"information is available"));
Assert (fe_index == active_fe_indices[obj_index],
- ExcMessage ("FE index does not match that of the present cell"));
+ ExcMessage ("FE index does not match that of the present cell"));
// see if the dof_indices array has been compressed for this
// particular cell
Assert ((signed_active_fe_index_type)active_fe_indices[obj_index]>=0,
ExcMessage ("This function can no longer be called after compressing the dof_indices array"));
Assert (fe_index == active_fe_indices[obj_index],
- ExcMessage ("FE index does not match that of the present cell"));
+ ExcMessage ("FE index does not match that of the present cell"));
dof_indices[dof_offsets[obj_index]+local_index] = global_index;
}
void
DoFLevel::
set_active_fe_index (const unsigned int obj_index,
- const unsigned int fe_index)
+ const unsigned int fe_index)
{
Assert (obj_index < active_fe_indices.size(),
ExcIndexRange (obj_index, 0, active_fe_indices.size()));
inline
const types::global_dof_index *
DoFLevel::get_cell_cache_start (const unsigned int obj_index,
- const unsigned int dofs_per_cell) const
+ const unsigned int dofs_per_cell) const
{
Assert (obj_index < cell_cache_offsets.size(),
ExcInternalError());
}
-/**
- * Cell matrix for gradient. The derivative is on the trial function.
- * \f[
- * \int_Z \nabla u \cdot \mathbf v\,dx
- * \f]
- *
- * This is the strong gradient and the trial space should be at least
- * in <i>H</i><sup>1</sup>. The test functions can be discontinuous.
- *
- * @author Guido Kanschat
- * @date 2011
- */
+ /**
+ * Cell matrix for gradient. The derivative is on the trial function.
+ * \f[
+ * \int_Z \nabla u \cdot \mathbf v\,dx
+ * \f]
+ *
+ * This is the strong gradient and the trial space should be at least
+ * in <i>H</i><sup>1</sup>. The test functions can be discontinuous.
+ *
+ * @author Guido Kanschat
+ * @date 2011
+ */
template <int dim>
void gradient_matrix(
FullMatrix<double> &M,
{
const double dx = fe1.JxW(k);
const Point<dim> &n = fe1.normal_vector(k);
- for (unsigned int i=0; i<n_dofs; ++i)
- {
- for (unsigned int j=0; j<n_dofs; ++j)
- {
- double u1dotn = 0.;
- double v1dotn = 0.;
- double u2dotn = 0.;
- double v2dotn = 0.;
-
- double ngradu1n = 0.;
- double ngradv1n = 0.;
- double ngradu2n = 0.;
- double ngradv2n = 0.;
-
- for (unsigned int d=0; d<dim; ++d)
- {
- u1dotn += n(d)*fe1.shape_value_component(j,k,d);
- v1dotn += n(d)*fe1.shape_value_component(i,k,d);
- u2dotn += n(d)*fe2.shape_value_component(j,k,d);
- v2dotn += n(d)*fe2.shape_value_component(i,k,d);
-
- ngradu1n += n*fe1.shape_grad_component(j,k,d)*n(d);
- ngradv1n += n*fe1.shape_grad_component(i,k,d)*n(d);
- ngradu2n += n*fe2.shape_grad_component(j,k,d)*n(d);
- ngradv2n += n*fe2.shape_grad_component(i,k,d)*n(d);
- }
-
- for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
- {
+ for (unsigned int i=0; i<n_dofs; ++i)
+ {
+ for (unsigned int j=0; j<n_dofs; ++j)
+ {
+ double u1dotn = 0.;
+ double v1dotn = 0.;
+ double u2dotn = 0.;
+ double v2dotn = 0.;
+
+ double ngradu1n = 0.;
+ double ngradv1n = 0.;
+ double ngradu2n = 0.;
+ double ngradv2n = 0.;
+
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ u1dotn += n(d)*fe1.shape_value_component(j,k,d);
+ v1dotn += n(d)*fe1.shape_value_component(i,k,d);
+ u2dotn += n(d)*fe2.shape_value_component(j,k,d);
+ v2dotn += n(d)*fe2.shape_value_component(i,k,d);
+
+ ngradu1n += n*fe1.shape_grad_component(j,k,d)*n(d);
+ ngradv1n += n*fe1.shape_grad_component(i,k,d)*n(d);
+ ngradu2n += n*fe2.shape_grad_component(j,k,d)*n(d);
+ ngradv2n += n*fe2.shape_grad_component(i,k,d)*n(d);
+ }
+
+ for (unsigned int d=0; d<fe1.get_fe().n_components(); ++d)
+ {
const double vi = fe1.shape_value_component(i,k,d)-v1dotn*n(d);
- const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n(d);
-
+ const double dnvi = n * fe1.shape_grad_component(i,k,d)-ngradv1n*n(d);
+
const double ve = fe2.shape_value_component(i,k,d)-v2dotn*n(d);
const double dnve = n * fe2.shape_grad_component(i,k,d)-ngradv2n*n(d);
-
+
const double ui = fe1.shape_value_component(j,k,d)-u1dotn*n(d);
const double dnui = n * fe1.shape_grad_component(j,k,d)-ngradu1n*n(d);
-
+
const double ue = fe2.shape_value_component(j,k,d)-u2dotn*n(d);
const double dnue = n * fe2.shape_grad_component(j,k,d)-ngradu2n*n(d);
- M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
+ M11(i,j) += dx*(-.5*nui*dnvi*ui-.5*nui*dnui*vi+nu*penalty*ui*vi);
M12(i,j) += dx*( .5*nui*dnvi*ue-.5*nue*dnue*vi-nu*penalty*vi*ue);
M21(i,j) += dx*(-.5*nue*dnve*ui+.5*nui*dnui*ve-nu*penalty*ui*ve);
M22(i,j) += dx*( .5*nue*dnve*ue+.5*nue*dnue*ve+nu*penalty*ue*ve);
* functions of ARPACK.
*/
template <typename VECTOR, typename MATRIX1,
- typename MATRIX2, typename INVERSE>
+ typename MATRIX2, typename INVERSE>
void solve(
const MATRIX1 &A,
const MATRIX2 &B,
template <typename VECTOR, typename MATRIX1,
- typename MATRIX2, typename INVERSE>
+ typename MATRIX2, typename INVERSE>
inline
void ArpackSolver::solve (
const MATRIX1 &system_matrix,
if (same_level_only && cell->neighbor(4)->neighbor(2)->level() != cell->level()) break;
if (cell->neighbor(4)->neighbor(0)->at_boundary(2)) break;
if (same_level_only && cell->neighbor(4)->neighbor(0)->neighbor(2)->level() != cell->level()) break;
- // No break here
+ // No break here
case 2:
if (cell->at_boundary(2)) break;
if (same_level_only && cell->neighbor(2)->level() != cell->level()) break;
indices[fe.face_to_cell_index(i,5)] = numbers::invalid_dof_index;
}
add(k, indices);
- // no break here
+ // no break here
case 1:
cell->get_mg_dof_indices(indices);
for (unsigned int i=0; i<fe.dofs_per_face; ++i)
* compiler.
*/
template <class BlockVectorType,
- class VectorType>
+ class VectorType>
void vmult_block_nonblock (BlockVectorType &dst,
const VectorType &src) const;
* compiler.
*/
template <class BlockVectorType,
- class VectorType>
+ class VectorType>
void vmult_nonblock_block (VectorType &dst,
const BlockVectorType &src) const;
* compiler.
*/
template <class BlockVectorType,
- class VectorType>
+ class VectorType>
void Tvmult_block_nonblock (BlockVectorType &dst,
const VectorType &src) const;
* compiler.
*/
template <class BlockVectorType,
- class VectorType>
+ class VectorType>
void Tvmult_nonblock_block (VectorType &dst,
const BlockVectorType &src) const;
* beginning of their use, so there is no point actually copying
* anything.
*/
- TemporaryData & operator = (const TemporaryData &)
+ TemporaryData &operator = (const TemporaryData &)
{
return *this;
}
template <class MatrixType>
template <class BlockVectorType,
- class VectorType>
+ class VectorType>
void
BlockMatrixBase<MatrixType>::
vmult_nonblock_block (VectorType &dst,
template <class MatrixType>
template <class BlockVectorType,
- class VectorType>
+ class VectorType>
void
BlockMatrixBase<MatrixType>::
vmult_block_nonblock (BlockVectorType &dst,
template <class MatrixType>
template <class BlockVectorType,
- class VectorType>
+ class VectorType>
void
BlockMatrixBase<MatrixType>::
Tvmult_block_nonblock (BlockVectorType &dst,
template <class MatrixType>
template <class BlockVectorType,
- class VectorType>
+ class VectorType>
void
BlockMatrixBase<MatrixType>::
Tvmult_nonblock_block (VectorType &dst,
* only one block column.
*/
template <typename block_number,
- typename nonblock_number>
+ typename nonblock_number>
void vmult (BlockVector<block_number> &dst,
const Vector<nonblock_number> &src) const;
* only one block row.
*/
template <typename block_number,
- typename nonblock_number>
+ typename nonblock_number>
void vmult (Vector<nonblock_number> &dst,
const BlockVector<block_number> &src) const;
* only one block row.
*/
template <typename block_number,
- typename nonblock_number>
+ typename nonblock_number>
void Tvmult (BlockVector<block_number> &dst,
const Vector<nonblock_number> &src) const;
* only one block column.
*/
template <typename block_number,
- typename nonblock_number>
+ typename nonblock_number>
void Tvmult (Vector<nonblock_number> &dst,
const BlockVector<block_number> &src) const;
template <typename number>
template <typename block_number,
- typename nonblock_number>
+ typename nonblock_number>
inline
void
BlockSparseMatrix<number>::vmult (BlockVector<block_number> &dst,
template <typename number>
template <typename block_number,
- typename nonblock_number>
+ typename nonblock_number>
inline
void
BlockSparseMatrix<number>::vmult (Vector<nonblock_number> &dst,
template <typename number>
template <typename block_number,
- typename nonblock_number>
+ typename nonblock_number>
inline
void
BlockSparseMatrix<number>::Tvmult (BlockVector<block_number> &dst,
template <typename number>
template <typename block_number,
- typename nonblock_number>
+ typename nonblock_number>
inline
void
BlockSparseMatrix<number>::Tvmult (Vector<nonblock_number> &dst,
template <typename OtherNumber>
inline
void BlockVectorBase<VectorType>::extract_subvector_to (const std::vector<size_type> &indices,
- std::vector<OtherNumber> &values) const
+ std::vector<OtherNumber> &values) const
{
for (size_type i = 0; i < indices.size(); ++i)
values[i] = operator()(indices[i]);
const ForwardIterator indices_end,
OutputIterator values_begin) const
{
- while (indices_begin != indices_end) {
- *values_begin = operator()(*indices_begin);
- indices_begin++; values_begin++;
- }
+ while (indices_begin != indices_end)
+ {
+ *values_begin = operator()(*indices_begin);
+ indices_begin++;
+ values_begin++;
+ }
}
#endif // DOXYGEN
* fragment.
*/
template <typename MatrixIterator,
- typename SrcIterator,
- typename DstIterator>
+ typename SrcIterator,
+ typename DstIterator>
inline
void
chunk_vmult_add (const size_type chunk_size,
* the residual.
*/
template <typename MatrixIterator,
- typename SrcIterator,
- typename DstIterator>
+ typename SrcIterator,
+ typename DstIterator>
inline
void
chunk_vmult_subtract (const size_type chunk_size,
* chunk_size to the destination vector fragment.
*/
template <typename MatrixIterator,
- typename SrcIterator,
- typename DstIterator>
+ typename SrcIterator,
+ typename DstIterator>
inline
void
chunk_Tvmult_add (const size_type chunk_size,
* individual chunk.
*/
template <typename result_type,
- typename MatrixIterator,
- typename SrcIterator1,
- typename SrcIterator2>
+ typename MatrixIterator,
+ typename SrcIterator1,
+ typename SrcIterator2>
inline
result_type
chunk_matrix_scalar_product (const size_type chunk_size,
* task scheduler.
*/
template <typename number,
- typename InVector,
- typename OutVector>
+ typename InVector,
+ typename OutVector>
void vmult_add_on_subrange (const ChunkSparsityPattern &cols,
const unsigned int begin_row,
const unsigned int end_row,
* IndexSet, whereas otherwise it contains the global problem size and the
* local range.
*/
- const IndexSet & get_local_lines() const;
+ const IndexSet &get_local_lines() const;
/**
* This function copies the content of @p constraints_in with DoFs that are
* <code>index2=index1</code>.
*/
bool are_identity_constrained (const size_type index1,
- const size_type index2) const;
+ const size_type index2) const;
/**
* Return the maximum number of other dofs that one dof is constrained
* prevent data races.
*/
template <typename ForwardIteratorVec, typename ForwardIteratorInd,
- class VectorType>
+ class VectorType>
void
distribute_local_to_global (ForwardIteratorVec local_vector_begin,
ForwardIteratorVec local_vector_end,
* set (i.e., distribute to be called).
*/
template <typename ForwardIteratorVec, typename ForwardIteratorInd,
- class VectorType>
+ class VectorType>
void
get_dof_values (const VectorType &global_vector,
ForwardIteratorInd local_indices_begin,
template <typename ForwardIteratorVec, typename ForwardIteratorInd,
- class VectorType>
+ class VectorType>
inline
void ConstraintMatrix::distribute_local_to_global (
ForwardIteratorVec local_vector_begin,
template <typename ForwardIteratorVec, typename ForwardIteratorInd,
- class VectorType>
+ class VectorType>
inline
void ConstraintMatrix::get_dof_values (const VectorType &global_vector,
ForwardIteratorInd local_indices_begin,
/**
* Dereferencing operator.
*/
- ScratchData & operator* ()
+ ScratchData &operator* ()
{
return *my_scratch_data;
}
/**
* Dereferencing operator.
*/
- ScratchData * operator-> ()
+ ScratchData *operator-> ()
{
return my_scratch_data;
}
private:
- ScratchData* my_scratch_data;
+ ScratchData *my_scratch_data;
};
private:
static Threads::ThreadLocalStorage<ScratchData> scratch_data;
};
-
+
// function for block matrices: Find out where in the list of local dofs
// (sorted according to global ids) the individual blocks start. Transform
const size_type n_local_dofs = local_dof_indices.size();
typename internals::ConstraintMatrixData<number>::ScratchDataAccessor
- scratch_data;
+ scratch_data;
internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
global_rows.reinit(n_local_dofs);
// an array in any case since we cannot know about the actual data type in
// the ConstraintMatrix class (unless we do cast). This involves a little
// bit of logic to determine the type of the matrix value.
- std::vector<size_type> & cols = scratch_data->columns;
- std::vector<number> & vals = scratch_data->values;
+ std::vector<size_type> &cols = scratch_data->columns;
+ std::vector<number> &vals = scratch_data->values;
SparseMatrix<number> *sparse_matrix
= dynamic_cast<SparseMatrix<number> *>(&global_matrix);
if (use_dealii_matrix == false)
const size_type n_local_col_dofs = col_indices.size();
typename internals::ConstraintMatrixData<number>::ScratchDataAccessor
- scratch_data;
+ scratch_data;
internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
global_rows.reinit(n_local_row_dofs);
internals::GlobalRowsFromLocal &global_cols = scratch_data->global_columns;
// create arrays for the column data (indices and values) that will then be
// written into the matrix. Shortcut for deal.II sparse matrix
- std::vector<size_type> & cols = scratch_data->columns;
- std::vector<number> & vals = scratch_data->values;
+ std::vector<size_type> &cols = scratch_data->columns;
+ std::vector<number> &vals = scratch_data->values;
cols.resize(n_actual_col_dofs);
vals.resize(n_actual_col_dofs);
Assert (sorted == true, ExcMatrixNotClosed());
typename internals::ConstraintMatrixData<number>::ScratchDataAccessor
- scratch_data;
+ scratch_data;
const size_type n_local_dofs = local_dof_indices.size();
internals::GlobalRowsFromLocal &global_rows = scratch_data->global_rows;
block_starts.resize(num_blocks+1);
internals::make_block_starts (global_matrix, global_rows, block_starts);
- std::vector<size_type> & cols = scratch_data->columns;
- std::vector<number> & vals = scratch_data->values;
+ std::vector<size_type> &cols = scratch_data->columns;
+ std::vector<number> &vals = scratch_data->values;
if (use_dealii_matrix == false)
{
cols.resize (n_actual_dofs);
// plus some indices that come from constraints.
if (dof_mask_is_active == false)
{
- std::vector<size_type> & actual_dof_indices = scratch_data->columns;
+ std::vector<size_type> &actual_dof_indices = scratch_data->columns;
actual_dof_indices.resize(n_local_dofs);
make_sorted_row_list (local_dof_indices, actual_dof_indices);
const size_type n_actual_dofs = actual_dof_indices.size();
else
sparsity_pattern.add (local_dof_indices[i], local_dof_indices[i]);
}
-
- return;
+
+ return;
}
// create arrays for the column indices that will then be written into the
// sparsity pattern.
- std::vector<size_type> & cols = scratch_data->columns;
+ std::vector<size_type> &cols = scratch_data->columns;
cols.resize(n_actual_dofs);
for (size_type i=0; i<n_actual_dofs; ++i)
if (dof_mask_is_active == false)
{
- std::vector<size_type> & actual_dof_indices = scratch_data->columns;
+ std::vector<size_type> &actual_dof_indices = scratch_data->columns;
actual_dof_indices.resize(n_local_dofs);
make_sorted_row_list (local_dof_indices, actual_dof_indices);
const size_type n_actual_dofs = actual_dof_indices.size();
const size_type n_actual_dofs = global_rows.size();
// additional construct that also takes care of block indices.
- std::vector<size_type> & block_starts = scratch_data->block_starts;
+ std::vector<size_type> &block_starts = scratch_data->block_starts;
block_starts.resize(num_blocks+1);
internals::make_block_starts(sparsity_pattern, global_rows, block_starts);
inline void
syev (const char *jobz, const char *uplo, const int *n, double *A, const int *lda, double *w, double *work, const int *lwork, int *info)
{
- dsyev_ (char*jobz,char*uplo,int*n,double*A,int*lda,double*w,double*work,int*lwork,int*info);
+ dsyev_ (char *jobz,char *uplo,int *n,double*A,int *lda,double*w,double*work,int *lwork,int *info);
}
#else
inline void
inline void
syev (const char *jobz, const char *uplo, const int *n, float *A, const int *lda, float *w, float *work, const int *lwork, int *info)
{
- ssyev_ (char*jobz,char*uplo,int*n,double*A,int*lda,double*w,double*work,int*lwork,int*info);
+ ssyev_ (char *jobz,char *uplo,int *n,double*A,int *lda,double*w,double*work,int *lwork,int *info);
}
#else
inline void
size_type n_elements = 0;
for (size_type row=i*options.block_size;
row < std::min(size_type(matrix.m()),
- size_type((i+1)*options.block_size)); ++row)
+ size_type((i+1)*options.block_size)); ++row)
for (size_type col=j*options.block_size;
col < std::min(size_type(matrix.m()),
- size_type((j+1)*options.block_size)); ++col, ++n_elements)
+ size_type((j+1)*options.block_size)); ++col, ++n_elements)
if (options.show_absolute_values == true)
average += std::fabs(get_element (matrix, row, col));
else
const Options options)
{
size_type
- gridpoints_x = (matrix.n() / options.block_size
- +
- (matrix.n() % options.block_size != 0 ? 1 : 0)),
- gridpoints_y = (matrix.m() / options.block_size
- +
- (matrix.m() % options.block_size != 0 ? 1 : 0));
+ gridpoints_x = (matrix.n() / options.block_size
+ +
+ (matrix.n() % options.block_size != 0 ? 1 : 0)),
+ gridpoints_y = (matrix.m() / options.block_size
+ +
+ (matrix.m() % options.block_size != 0 ? 1 : 0));
// If continuous, the number of
// plotted patches is matrix size-1
template <typename Number>
inline
BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
- const std::vector<IndexSet> &ghost_indices,
- const MPI_Comm communicator)
+ const std::vector<IndexSet> &ghost_indices,
+ const MPI_Comm communicator)
{
std::vector<size_type> sizes(local_ranges.size());
for (unsigned int i=0; i<local_ranges.size(); ++i)
template <typename Number>
inline
BlockVector<Number>::BlockVector (const std::vector<IndexSet> &local_ranges,
- const MPI_Comm communicator)
+ const MPI_Comm communicator)
{
std::vector<size_type> sizes(local_ranges.size());
for (unsigned int i=0; i<local_ranges.size(); ++i)
// we only allow assignment to vectors with the same number of blocks
// or to an empty BlockVector
Assert (this->n_blocks() == 0 || this->n_blocks() == v.n_blocks(),
- ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
+ ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
if (this->n_blocks() != v.n_blocks())
reinit(v.n_blocks(), true);
while (indices_begin != indices_end)
{
*values_begin = operator()(*indices_begin);
- indices_begin++; values_begin++;
+ indices_begin++;
+ values_begin++;
}
}
const size_type row,
const size_type index);
- /**
- * Copy constructor.
- */
- Accessor (const Accessor &a);
+ /**
+ * Copy constructor.
+ */
+ Accessor (const Accessor &a);
/**
* Row number of the element
// we only allow assignment to vectors with the same number of blocks
// or to an empty BlockVector
Assert (n_blocks() == 0 || n_blocks() == v.n_blocks(),
- ExcDimensionMismatch(n_blocks(), v.n_blocks()));
+ ExcDimensionMismatch(n_blocks(), v.n_blocks()));
if (this->n_blocks() != v.n_blocks())
reinit(v.n_blocks());
// if we are dealing
// with a parallel vector
if (ghosted )
- {
-
- int ierr;
-
- // there is the possibility
- // that the vector has
- // ghost elements. in that
- // case, we first need to
- // figure out which
- // elements we own locally,
- // then get a pointer to
- // the elements that are
- // stored here (both the
- // ones we own as well as
- // the ghost elements). in
- // this array, the locally
- // owned elements come
- // first followed by the
- // ghost elements whose
- // position we can get from
- // an index set
- PetscInt begin, end, i;
- ierr = VecGetOwnershipRange (vector, &begin, &end);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
+ {
+
+ int ierr;
+
+ // there is the possibility
+ // that the vector has
+ // ghost elements. in that
+ // case, we first need to
+ // figure out which
+ // elements we own locally,
+ // then get a pointer to
+ // the elements that are
+ // stored here (both the
+ // ones we own as well as
+ // the ghost elements). in
+ // this array, the locally
+ // owned elements come
+ // first followed by the
+ // ghost elements whose
+ // position we can get from
+ // an index set
+ PetscInt begin, end, i;
+ ierr = VecGetOwnershipRange (vector, &begin, &end);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ Vec locally_stored_elements = PETSC_NULL;
+ ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ PetscInt lsize;
+ ierr = VecGetSize(locally_stored_elements, &lsize);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ PetscScalar *ptr;
+ ierr = VecGetArray(locally_stored_elements, &ptr);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ for (i = 0; i < n_idx; i++)
+ {
+ const unsigned int index = *(indices_begin+i);
+ if ( index>=static_cast<unsigned int>(begin)
+ && index<static_cast<unsigned int>(end) )
+ {
+ //local entry
+ *(values_begin+i) = *(ptr+index-begin);
+ }
+ else
+ {
+ //ghost entry
+ const unsigned int ghostidx
+ = ghost_indices.index_within_set(index);
+
+ Assert(ghostidx+end-begin<(unsigned int)lsize, ExcInternalError());
+ *(values_begin+i) = *(ptr+ghostidx+end-begin);
+ }
+ }
+
+ ierr = VecRestoreArray(locally_stored_elements, &ptr);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
+
+ ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
- Vec locally_stored_elements = PETSC_NULL;
- ierr = VecGhostGetLocalForm(vector, &locally_stored_elements);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- PetscInt lsize;
- ierr = VecGetSize(locally_stored_elements, &lsize);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- PetscScalar *ptr;
- ierr = VecGetArray(locally_stored_elements, &ptr);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- for (i = 0; i < n_idx; i++) {
- const unsigned int index = *(indices_begin+i);
- if ( index>=static_cast<unsigned int>(begin)
- && index<static_cast<unsigned int>(end) )
- {
- //local entry
- *(values_begin+i) = *(ptr+index-begin);
- }
- else
- {
- //ghost entry
- const unsigned int ghostidx
- = ghost_indices.index_within_set(index);
-
- Assert(ghostidx+end-begin<(unsigned int)lsize, ExcInternalError());
- *(values_begin+i) = *(ptr+ghostidx+end-begin);
- }
}
-
- ierr = VecRestoreArray(locally_stored_elements, &ptr);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- ierr = VecGhostRestoreLocalForm(vector, &locally_stored_elements);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
-
- }
// if the vector is local or the
// caller, then simply access the
// element we are interested in
else
- {
- int ierr;
+ {
+ int ierr;
- PetscInt begin, end;
- ierr = VecGetOwnershipRange (vector, &begin, &end);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
+ PetscInt begin, end;
+ ierr = VecGetOwnershipRange (vector, &begin, &end);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
- PetscScalar *ptr;
- ierr = VecGetArray(vector, &ptr);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
+ PetscScalar *ptr;
+ ierr = VecGetArray(vector, &ptr);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
- for (PetscInt i = 0; i < n_idx; i++) {
- const unsigned int index = *(indices_begin+i);
+ for (PetscInt i = 0; i < n_idx; i++)
+ {
+ const unsigned int index = *(indices_begin+i);
- Assert(index>=static_cast<unsigned int>(begin)
- && index<static_cast<unsigned int>(end), ExcInternalError());
+ Assert(index>=static_cast<unsigned int>(begin)
+ && index<static_cast<unsigned int>(end), ExcInternalError());
- *(values_begin+i) = *(ptr+index-begin);
- }
+ *(values_begin+i) = *(ptr+index-begin);
+ }
- ierr = VecRestoreArray(vector, &ptr);
- AssertThrow (ierr == 0, ExcPETScError(ierr));
+ ierr = VecRestoreArray(vector, &ptr);
+ AssertThrow (ierr == 0, ExcPETScError(ierr));
- }
+ }
}
#endif // DOXYGEN
* @author Ralf Hartmann, 1999
*/
template <class MATRIX = SparseMatrix<double>,
- class VECTOR = dealii::Vector<double> >
+ class VECTOR = dealii::Vector<double> >
class PreconditionSelector : public Subscriptor
{
public:
*/
void
solve (const unsigned int n_eigenpairs,
- unsigned int *n_converged);
+ unsigned int *n_converged);
/**
* Access the real parts of solutions for a solved eigenvector
* preconditioners.
*/
class AdditionalData
-{};
+ {};
/**
*/
template <class STREAM>
void print (STREAM &out,
- const bool across = false,
- const bool diagonal_first = true) const;
+ const bool across = false,
+ const bool diagonal_first = true) const;
/**
* Print the matrix in the usual format, i.e. as a matrix and not as a list
template <class STREAM>
inline
void SparseMatrix<number>::print (STREAM &out,
- const bool across,
- const bool diagonal_first) const
+ const bool across,
+ const bool diagonal_first) const
{
Assert (cols != 0, ExcNotInitialized());
Assert (val != 0, ExcNotInitialized());
* task scheduler.
*/
template <typename number,
- typename InVector,
- typename OutVector>
+ typename InVector,
+ typename OutVector>
void vmult_on_subrange (const size_type begin_row,
const size_type end_row,
const number *values,
* scheduler.
*/
template <typename number,
- typename InVector>
+ typename InVector>
number matrix_norm_sqr_on_subrange (const size_type begin_row,
const size_type end_row,
const number *values,
* scheduler.
*/
template <typename number,
- typename InVector>
+ typename InVector>
number matrix_scalar_product_on_subrange (const size_type begin_row,
const size_type end_row,
const number *values,
* scheduler.
*/
template <typename number,
- typename InVector,
- typename OutVector>
+ typename InVector,
+ typename OutVector>
number residual_sqr_on_subrange (const size_type begin_row,
const size_type end_row,
const number *values,
/**
* similar to the function above, but includes support for
* BlockCompressedSimpleSparsityPattern.
- * @p owned_set_per_cpu is typically DoFHandler::locally_owned_dofs_per_processor
+ * @p owned_set_per_cpu is typically DoFHandler::locally_owned_dofs_per_processor
* and @p myrange are locally_relevant_dofs.
*/
template <class CSP_t>
* be ignored.
*/
BlockVector (const std::vector<IndexSet> ¶llel_partitioning,
- const std::vector<IndexSet> &ghost_values,
- const MPI_Comm &communicator);
+ const std::vector<IndexSet> &ghost_values,
+ const MPI_Comm &communicator);
/**
inline
BlockVector::BlockVector (const std::vector<IndexSet> ¶llel_partitioning,
- const std::vector<IndexSet> &ghost_values,
- const MPI_Comm &communicator)
+ const std::vector<IndexSet> &ghost_values,
+ const MPI_Comm &communicator)
{
reinit(parallel_partitioning, ghost_values, communicator);
}
* Apply the transpose preconditioner.
*/
virtual void Tvmult (VectorBase &dst,
- const VectorBase &src) const;
+ const VectorBase &src) const;
/**
* Apply the preconditioner on
const VectorBase &src) const;
/**
- * Apply the preconditioner on deal.II data structures
- * instead of the ones provided in the Trilinos wrapper class,
+ * Apply the preconditioner on deal.II data structures
+ * instead of the ones provided in the Trilinos wrapper class,
* i.e., dst = src.
*/
void vmult (dealii::Vector<double> &dst,
const dealii::Vector<double> &src) const;
/**
- * Apply the transpose preconditioner on deal.II data structures
+ * Apply the transpose preconditioner on deal.II data structures
* instead of the ones provided in the Trilinos wrapper class,
* i.e. dst = src.
*/
const dealii::Vector<double> &src) const;
/**
- * Apply the preconditioner on deal.II parallel data structures
+ * Apply the preconditioner on deal.II parallel data structures
* instead of the ones provided in the Trilinos wrapper class,
* i.e., dst = src.
*/
const dealii::parallel::distributed::Vector<double> &src) const;
/**
- * Apply the transpose preconditioner on deal.II parallel data structures
+ * Apply the transpose preconditioner on deal.II parallel data structures
* instead of the ones provided in the Trilinos wrapper class,
* i.e., dst = src.
*/
const ForwardIterator indices_end,
OutputIterator values_begin) const
{
- while (indices_begin != indices_end) {
- *values_begin = operator()(*indices_begin);
- indices_begin++; values_begin++;
- }
+ while (indices_begin != indices_end)
+ {
+ *values_begin = operator()(*indices_begin);
+ indices_begin++;
+ values_begin++;
+ }
}
}
-
+
#endif // DOXYGEN
}
const ForwardIterator indices_end,
OutputIterator values_begin) const
{
- while (indices_begin != indices_end) {
- *values_begin = operator()(*indices_begin);
- indices_begin++; values_begin++;
- }
+ while (indices_begin != indices_end)
+ {
+ *values_begin = operator()(*indices_begin);
+ indices_begin++;
+ values_begin++;
+ }
}
// The code returns the result as the last argument in order to make
// spawning tasks simpler and use automatic template deduction.
template <typename Operation, typename Number, typename Number2,
- typename ResultType, typename size_type>
+ typename ResultType, typename size_type>
void accumulate (const Operation &op,
const Number *X,
const Number2 *Y,
r2 = op(X, Y, power);
for (size_type j=1; j<8; ++j)
r2 += op(X, Y, power);
- // no break
+ // no break
case 2:
r1 = op(X, Y, power);
for (size_type j=1; j<8; ++j)
r1 += op(X, Y, power);
r1 += r2;
- // no break
+ // no break
case 1:
r2 = op(X, Y, power);
for (size_type j=1; j<8; ++j)
r2 += op(X, Y, power);
- // no break
+ // no break
default:
for (size_type j=0; j<remainder_inner; ++j)
r0 += op(X, Y, power);
// locally owned plus the number of ghosts
const unsigned int index_range = (vector_partitioner->local_range().second-
vector_partitioner->local_range().first)
- + vector_partitioner->ghost_indices().n_elements();
+ + vector_partitioner->ghost_indices().n_elements();
for (std::size_t i=0; i<dof_indices.size(); ++i)
AssertIndexRange (dof_indices[i], index_range);
else
{
const unsigned int *it = begin_indices (block),
- * end_cell = end_indices (block);
+ * end_cell = end_indices (block);
for ( ; it != end_cell; ++it)
if (row_lengths[*it] > 0)
{
{
out << "Entries row " << row << ": ";
const unsigned int *glob_indices = begin_indices(row),
- *end_row = end_indices(row);
+ *end_row = end_indices(row);
unsigned int index = 0;
const std::pair<unsigned short,unsigned short>
*con_it = begin_indicators(row),
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
class FEEvaluationBase
{
public:
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
class FEEvaluationAccess :
public FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
{
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
- int n_components_ = 1, typename Number = double >
+ int n_components_ = 1, typename Number = double >
class FEEvaluationGeneral :
public FEEvaluationAccess<dim,
Utilities::fixed_int_power<fe_degree+1,dim>::value,
* @author Katharina Kormann and Martin Kronbichler, 2010, 2011
*/
template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
- int n_components_ = 1, typename Number = double >
+ int n_components_ = 1, typename Number = double >
class FEEvaluation :
public FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
{
/*----------------------- FEEvaluationBase ----------------------------------*/
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
::FEEvaluationBase (const MatrixFree<dim,Number> &data_in,
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
void
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
unsigned int
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
internal::MatrixFreeFunctions::CellType
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
VectorType &vec,
Number &res) const
{
- res = vector_access (const_cast<const VectorType&>(vec), index);
+ res = vector_access (const_cast<const VectorType &>(vec), index);
}
void pre_constraints (const Number &,
VectorType &vec,
Number &res) const
{
- res += weight * vector_access (const_cast<const VectorType&>(vec), index);
+ res += weight * vector_access (const_cast<const VectorType &>(vec), index);
}
void post_constraints (const Number &sum,
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType, typename VectorOperation>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
template<typename VectorType>
inline
void
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components, typename Number>
inline
const VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components, typename Number>
inline
VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components, typename Number>
inline
const VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components, typename Number>
inline
VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components, typename Number>
inline
const VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components, typename Number>
inline
VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components, typename Number>
inline
const VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components, typename Number>
+ int n_components, typename Number>
inline
VectorizedArray<Number> *
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components,Number>::
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
Tensor<1,n_components_,VectorizedArray<Number> >
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
Tensor<1,n_components_,VectorizedArray<Number> >
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > >
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
Tensor<1,n_components_,Tensor<2,dim,VectorizedArray<Number> > >
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > >
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
Tensor<1,n_components_,VectorizedArray<Number> >
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
void
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
void
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
void
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
Tensor<1,n_components_,VectorizedArray<Number> >
FEEvaluationBase<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
template <int dim, int dofs_per_cell_, int n_q_points_,
- int n_components_, typename Number>
+ int n_components_, typename Number>
inline
FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,n_components_,Number>
::FEEvaluationAccess (const MatrixFree<dim,Number> &data_in,
}
template <int dim, int dofs_per_cell_, int n_q_points_,
- typename Number>
+ typename Number>
inline
void
FEEvaluationAccess<dim,dofs_per_cell_,n_q_points_,dim,Number>
/*----------------------- FEEvaluationGeneral -------------------------------*/
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
inline
FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
::FEEvaluationGeneral (const MatrixFree<dim,Number> &data_in,
{
proposed_quad_comp = no;
break;
- }
+ }
if (proposed_dof_comp != numbers::invalid_unsigned_int &&
proposed_quad_comp != numbers::invalid_unsigned_int)
{
// like the functions below and corresponds to a usual matrix-matrix
// product
template <int dim, int fe_degree, int n_q_points_1d, typename Number,
- int direction, bool dof_to_quad, bool add>
+ int direction, bool dof_to_quad, bool add>
inline
void
apply_tensor_product (const Number *shape_data,
// fe_degree+1 degrees of freedom per direction and not n_q_points_1d for
// those directions lower than the one currently applied
template <int dim, int fe_degree, typename Number, int face_direction,
- bool dof_to_quad, bool add>
+ bool dof_to_quad, bool add>
inline
void
apply_tensor_product_face (const Number *shape_data,
// ones and in addition use the symmetry in entries to reduce the number of
// read operations.
template <int dim, int fe_degree, int n_q_points_1d, typename Number,
- int direction, bool dof_to_quad, bool add>
+ int direction, bool dof_to_quad, bool add>
inline
void
apply_tensor_product_values (const Number *shape_values,
// zeros and ones and in addition use the symmetry in entries to
// reduce the number of read operations.
template <int dim, int fe_degree, int n_q_points_1d, typename Number,
- int direction, bool dof_to_quad, bool add>
+ int direction, bool dof_to_quad, bool add>
inline
void
apply_tensor_product_gradients (const Number *shape_gradients,
// form assuming the symmetries of unit cell shape hessians for
// finite elements in FEEvaluation
template <int dim, int fe_degree, int n_q_points_1d, typename Number,
- int direction, bool dof_to_quad, bool add>
+ int direction, bool dof_to_quad, bool add>
inline
void
apply_tensor_product_hessians (const Number *shape_hessians,
// [ 2.82 -1.528 2.673 0 -13.51 ]
// [-1 0.518 -0.75 2.482 10 ]
template <int dim, int fe_degree, typename Number,
- int direction, bool dof_to_quad, bool add>
+ int direction, bool dof_to_quad, bool add>
inline
void
apply_tensor_product_gradients_gl (const Number *shape_gradients,
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
inline
void
FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
inline
void
FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
inline
Point<dim,VectorizedArray<Number> >
FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
inline
FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
::FEEvaluation (const MatrixFree<dim,Number> &data_in,
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
inline
void
FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
inline
void
FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
// for all other orders, we choose even-odd
if (fe_degree > 1 || n_q_points_1d > 3)
internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add, 0>
- (shape_val_evenodd, in, out);
+ VectorizedArray<Number>, direction, dof_to_quad, add, 0>
+ (shape_val_evenodd, in, out);
else
internal::apply_tensor_product_values<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add>
- (this->data.shape_values.begin(), in, out);
+ VectorizedArray<Number>, direction, dof_to_quad, add>
+ (this->data.shape_values.begin(), in, out);
}
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
{
if (fe_degree > 1 || n_q_points_1d > 3)
internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add, 1>
- (shape_gra_evenodd, in, out);
+ VectorizedArray<Number>, direction, dof_to_quad, add, 1>
+ (shape_gra_evenodd, in, out);
else
internal::apply_tensor_product_gradients<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add>
- (this->data.shape_gradients.begin(), in, out);
+ VectorizedArray<Number>, direction, dof_to_quad, add>
+ (this->data.shape_gradients.begin(), in, out);
}
// the same symmetry relations hold. However, it is not possible to omit some
// values that are zero for the values
template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
- typename Number>
+ typename Number>
template <int direction, bool dof_to_quad, bool add>
inline
void
{
if (fe_degree > 1 || n_q_points_1d > 3)
internal::apply_tensor_product_evenodd<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add, 2>
- (shape_hes_evenodd, in, out);
+ VectorizedArray<Number>, direction, dof_to_quad, add, 2>
+ (shape_hes_evenodd, in, out);
else
internal::apply_tensor_product_hessians<dim,fe_degree,n_q_points_1d,
- VectorizedArray<Number>, direction, dof_to_quad, add>
- (this->data.shape_hessians.begin(), in, out);
+ VectorizedArray<Number>, direction, dof_to_quad, add>
+ (this->data.shape_hessians.begin(), in, out);
}
for (types::global_dof_index i=start_index; i<end_index; ++i)
if (constraint[no]->is_constrained(i)==true)
dof_info[no].constrained_dofs.
- push_back(static_cast<unsigned int>(i-start_index));
+ push_back(static_cast<unsigned int>(i-start_index));
}
}
const std_cxx1x::function<void (DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker,
const std_cxx1x::function<void (DOFINFO &, DOFINFO &,
typename INFOBOX::CellInfo &,
- typename INFOBOX::CellInfo &)>& face_worker,
+ typename INFOBOX::CellInfo &)> &face_worker,
const bool cells_first,
const bool unique_faces_only)
{
const std_cxx1x::function<void (DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker,
const std_cxx1x::function<void (DOFINFO &, DOFINFO &,
typename INFOBOX::CellInfo &,
- typename INFOBOX::CellInfo &)>& face_worker,
+ typename INFOBOX::CellInfo &)> &face_worker,
ASSEMBLER &assembler,
bool cells_first = true,
bool unique_faces_only = true)
for (unsigned int j=0; j<i1.size(); ++j)
for (unsigned int k=0; k<i2.size(); ++k)
if (std::fabs(M(k,j)) >= threshold)
- if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+ if (!mg_constrained_dofs->at_refinement_edge(level, i2[k]))
G.add(i1[j], i2[k], M(k,j));
}
}
AssertDimension(M.m(), i1.size());
AssertDimension(M.n(), i2.size());
Assert(mg_constrained_dofs != 0, ExcInternalError());
-
+
for (unsigned int j=0; j<i1.size(); ++j)
- for (unsigned int k=0; k<i2.size(); ++k)
- if (std::fabs(M(j,k)) >= threshold)
- // Enter values into matrix only if j corresponds to a
- // degree of freedom on the refinemenent edge, k does
- // not, and both are not on the boundary. This is part
- // the difference between the complete matrix with no
- // boundary condition at the refinement edge and and
- // the matrix assembled above by assemble().
-
- // Thus the logic is: enter the row if it is
- // constrained by hanging node constraints (actually,
- // the whole refinement edge), but not if it is
- // constrained by a boundary constraint.
- if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
- !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
- {
- if (mg_constrained_dofs->set_boundary_values())
- {
- if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
- !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
- ||
- (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
- mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
- i1[j] == i2[k]))
- G.add(i1[j], i2[k], M(j,k));
- }
- else
- G.add(i1[j], i2[k], M(j,k));
- }
+ for (unsigned int k=0; k<i2.size(); ++k)
+ if (std::fabs(M(j,k)) >= threshold)
+ // Enter values into matrix only if j corresponds to a
+ // degree of freedom on the refinemenent edge, k does
+ // not, and both are not on the boundary. This is part
+ // the difference between the complete matrix with no
+ // boundary condition at the refinement edge and and
+ // the matrix assembled above by assemble().
+
+ // Thus the logic is: enter the row if it is
+ // constrained by hanging node constraints (actually,
+ // the whole refinement edge), but not if it is
+ // constrained by a boundary constraint.
+ if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+ {
+ if (mg_constrained_dofs->set_boundary_values())
+ {
+ if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+ ||
+ (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+ mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+ i1[j] == i2[k]))
+ G.add(i1[j], i2[k], M(j,k));
+ }
+ else
+ G.add(i1[j], i2[k], M(j,k));
+ }
}
-
-
+
+
template <class MATRIX>
inline void
MGMatrixSimple<MATRIX>::assemble_out(
AssertDimension(M.n(), i1.size());
AssertDimension(M.m(), i2.size());
Assert(mg_constrained_dofs != 0, ExcInternalError());
-
+
for (unsigned int j=0; j<i1.size(); ++j)
- for (unsigned int k=0; k<i2.size(); ++k)
- if (std::fabs(M(k,j)) >= threshold)
- if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
- !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
- {
- if (mg_constrained_dofs->set_boundary_values())
- {
- if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
- !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
- ||
- (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
- mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
- i1[j] == i2[k]))
- G.add(i1[j], i2[k], M(k,j));
- }
- else
- G.add(i1[j], i2[k], M(k,j));
- }
+ for (unsigned int k=0; k<i2.size(); ++k)
+ if (std::fabs(M(k,j)) >= threshold)
+ if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+ {
+ if (mg_constrained_dofs->set_boundary_values())
+ {
+ if ((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+ !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+ ||
+ (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+ mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+ i1[j] == i2[k]))
+ G.add(i1[j], i2[k], M(k,j));
+ }
+ else
+ G.add(i1[j], i2[k], M(k,j));
+ }
}
-
+
template <class MATRIX>
template <class DOFINFO>
ResidualSimple<VECTOR>::initialize(c);
}
-
+
template <class MATRIX, class VECTOR>
template <class DOFINFO>
inline void
namespace mg
{
-/**
- * Handler and storage for all five SparseMatrix object involved in
- * using multigrid with local refinement.
- *
- * @author Baerbel Janssen, Guido Kanschat
- * @date 2013
- */
+ /**
+ * Handler and storage for all five SparseMatrix object involved in
+ * using multigrid with local refinement.
+ *
+ * @author Baerbel Janssen, Guido Kanschat
+ * @date 2013
+ */
template <typename number>
class SparseMatrixCollection : public Subscriptor
{
- public:
- void resize(const unsigned int minlevel, const unsigned int maxlevel);
-
- template <class DH>
- void reinit(const DH& dof_handler);
-
- void set_zero();
-
- MGLevelObject<SparsityPattern> sparsity;
- MGLevelObject<SparsityPattern> sparsity_edge;
-
- MGLevelObject<SparseMatrix<number> > matrix;
- MGLevelObject<SparseMatrix<number> > matrix_down;
- MGLevelObject<SparseMatrix<number> > matrix_up;
- MGLevelObject<SparseMatrix<number> > matrix_in;
- MGLevelObject<SparseMatrix<number> > matrix_out;
+ public:
+ void resize(const unsigned int minlevel, const unsigned int maxlevel);
+
+ template <class DH>
+ void reinit(const DH &dof_handler);
+
+ void set_zero();
+
+ MGLevelObject<SparsityPattern> sparsity;
+ MGLevelObject<SparsityPattern> sparsity_edge;
+
+ MGLevelObject<SparseMatrix<number> > matrix;
+ MGLevelObject<SparseMatrix<number> > matrix_down;
+ MGLevelObject<SparseMatrix<number> > matrix_up;
+ MGLevelObject<SparseMatrix<number> > matrix_in;
+ MGLevelObject<SparseMatrix<number> > matrix_out;
};
template <typename number>
template <class DH>
void
- SparseMatrixCollection<number>::reinit(const DH& dof_handler)
+ SparseMatrixCollection<number>::reinit(const DH &dof_handler)
{
AssertIndexRange(sparsity.max_level(), dof_handler.get_tria().n_levels());
-
+
for (unsigned int level=sparsity.min_level();
- level<=sparsity.max_level();++level)
+ level<=sparsity.max_level(); ++level)
{
- CompressedSparsityPattern c_sparsity(dof_handler.n_dofs(level));
- MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level);
- sparsity[level].copy_from(c_sparsity);
- matrix[level].reinit(sparsity[level]);
- matrix_in[level].reinit(sparsity[level]);
- matrix_out[level].reinit(sparsity[level]);
- if (level>0)
- {
- CompressedSparsityPattern ci_sparsity;
- ci_sparsity.reinit(dof_handler.n_dofs(level-1), dof_handler.n_dofs(level));
- MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level);
- sparsity_edge[level].copy_from(ci_sparsity);
- matrix_up[level].reinit(sparsity_edge[level]);
- matrix_down[level].reinit(sparsity_edge[level]);
- }
+ CompressedSparsityPattern c_sparsity(dof_handler.n_dofs(level));
+ MGTools::make_flux_sparsity_pattern(dof_handler, c_sparsity, level);
+ sparsity[level].copy_from(c_sparsity);
+ matrix[level].reinit(sparsity[level]);
+ matrix_in[level].reinit(sparsity[level]);
+ matrix_out[level].reinit(sparsity[level]);
+ if (level>0)
+ {
+ CompressedSparsityPattern ci_sparsity;
+ ci_sparsity.reinit(dof_handler.n_dofs(level-1), dof_handler.n_dofs(level));
+ MGTools::make_flux_sparsity_pattern_edge(dof_handler, ci_sparsity, level);
+ sparsity_edge[level].copy_from(ci_sparsity);
+ matrix_up[level].reinit(sparsity_edge[level]);
+ matrix_down[level].reinit(sparsity_edge[level]);
+ }
}
}
matrix_up = 0.;
matrix_down = 0.;
}
-
+
}
DEAL_II_NAMESPACE_CLOSE
* we are to work on.
*/
template <class DerivativeDescription, int dim,
- template <int, int> class DH, class InputVector, int spacedim>
+ template <int, int> class DH, class InputVector, int spacedim>
static void
approximate_derivative (const Mapping<dim,spacedim> &mapping,
const DH<dim,spacedim> &dof,
* on the cell.
*/
template <class DerivativeDescription, int dim,
- template <int, int> class DH, class InputVector, int spacedim>
+ template <int, int> class DH, class InputVector, int spacedim>
static void
approximate (SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DH<dim,spacedim>::active_cell_iterator>,
- Vector<float>::iterator> > const &cell,
+ Vector<float>::iterator> > const &cell,
const Mapping<dim,spacedim> &mapping,
const DH<dim,spacedim> &dof,
const InputVector &solution,
* derivative tensor.
*/
template <class DerivativeDescription, int dim,
- template <int, int> class DH, class InputVector, int spacedim>
+ template <int, int> class DH, class InputVector, int spacedim>
static void
approximate_cell (const Mapping<dim,spacedim> &mapping,
const DH<dim,spacedim> &dof,
* @author Luca Heltai, 2006, Markus Buerg, 2012, Wolfgang Bangerth, 2013
*/
template <int dim,
- typename DH=DoFHandler<dim>,
- typename VECTOR=Vector<double> >
+ typename DH=DoFHandler<dim>,
+ typename VECTOR=Vector<double> >
class FEFieldFunction : public Function<dim>
{
public:
* constraints object.
*/
template <int dim, int spacedim,
- template <int,int> class DH,
- class VECTOR>
+ template <int,int> class DH,
+ class VECTOR>
void
interpolate_to_different_mesh (const DH<dim, spacedim> &dof1,
const VECTOR &u1,
* at hanging nodes.
*/
template <int dim, int spacedim,
- template <int,int> class DH,
- class VECTOR>
+ template <int,int> class DH,
+ class VECTOR>
void
interpolate_to_different_mesh (const DH<dim, spacedim> &dof1,
const VECTOR &u1,
* destination DoFHandler.
*/
template <int dim, int spacedim,
- template <int,int> class DH,
- class VECTOR>
+ template <int,int> class DH,
+ class VECTOR>
void
interpolate_to_different_mesh (const InterGridMap<DH<dim, spacedim> > &intergridmap,
const VECTOR &u1,
VECTOR &vec,
const bool enforce_zero_boundary = false,
const hp::QCollection<dim-1> &q_boundary = hp::QCollection<dim-1>(dim > 1 ?
- QGauss<dim-1>(2) :
- Quadrature<dim-1>(0)),
+ QGauss<dim-1>(2) :
+ Quadrature<dim-1>(0)),
const bool project_to_boundary_first = false);
/**
VECTOR &vec,
const bool enforce_zero_boundary = false,
const hp::QCollection<dim-1> &q_boundary = hp::QCollection<dim-1>(dim > 1 ?
- QGauss<dim-1>(2) :
- Quadrature<dim-1>(0)),
+ QGauss<dim-1>(2) :
+ Quadrature<dim-1>(0)),
const bool project_to_boundary_first = false);
/**
const bool fe_is_system = (n_components != 1);
typename DH<dim,spacedim>::active_cell_iterator cell = dof.begin_active(),
- endc = dof.end();
+ endc = dof.end();
// For FESystems many of the
// unit_support_points will appear
hp::MappingCollection<dim,spacedim> mapping_collection (mapping);
hp::FEValues<dim,spacedim> fe_values (mapping_collection,
- fe, support_quadrature, update_quadrature_points);
+ fe, support_quadrature, update_quadrature_points);
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
template <int dim, int spacedim,
- template <int,int> class DH,
- class Vector>
+ template <int,int> class DH,
+ class Vector>
void
interpolate_to_different_mesh (const DH<dim, spacedim> &dof1,
const Vector &u1,
template <int dim, int spacedim,
- template <int,int> class DH,
- class Vector>
+ template <int,int> class DH,
+ class Vector>
void
interpolate_to_different_mesh (const DH<dim, spacedim> &dof1,
const Vector &u1,
template <int dim, int spacedim,
- template <int,int> class DH,
- class Vector>
+ template <int,int> class DH,
+ class Vector>
void
interpolate_to_different_mesh (const InterGridMap<DH<dim, spacedim> > &intergridmap,
const Vector &u1,
template <int,int> class M_or_MC,
template <int> class Q_or_QC>
void project_compute_b_v (const M_or_MC<dim, spacedim> &mapping,
- const DH<dim,spacedim> &dof,
- const Function<spacedim> &function,
- const bool enforce_zero_boundary,
- const Q_or_QC<dim-1> &q_boundary,
- const bool project_to_boundary_first,
- std::map<types::global_dof_index,double> &boundary_values)
+ const DH<dim,spacedim> &dof,
+ const Function<spacedim> &function,
+ const bool enforce_zero_boundary,
+ const Q_or_QC<dim-1> &q_boundary,
+ const bool project_to_boundary_first,
+ std::map<types::global_dof_index,double> &boundary_values)
{
if (enforce_zero_boundary == true)
// no need to project boundary
template <int,int> class M_or_MC,
template <int> class Q_or_QC>
void do_project (const M_or_MC<dim, spacedim> &mapping,
- const DH<dim,spacedim> &dof,
- const ConstraintMatrix &constraints,
- const Q_or_QC<dim> &quadrature,
- const Function<spacedim> &function,
- Vector &vec_result,
- const bool enforce_zero_boundary,
- const Q_or_QC<dim-1> &q_boundary,
- const bool project_to_boundary_first)
+ const DH<dim,spacedim> &dof,
+ const ConstraintMatrix &constraints,
+ const Q_or_QC<dim> &quadrature,
+ const Function<spacedim> &function,
+ Vector &vec_result,
+ const bool enforce_zero_boundary,
+ const Q_or_QC<dim-1> &q_boundary,
+ const bool project_to_boundary_first)
{
Assert (dof.get_fe().n_components() == function.n_components,
ExcDimensionMismatch(dof.get_fe().n_components(),
const bool project_to_boundary_first)
{
do_project (mapping, dof, constraints, quadrature,
- function, vec_result,
- enforce_zero_boundary, q_boundary,
- project_to_boundary_first);
+ function, vec_result,
+ enforce_zero_boundary, q_boundary,
+ project_to_boundary_first);
}
const bool project_to_boundary_first)
{
do_project (mapping, dof, constraints, quadrature,
- function, vec_result,
- enforce_zero_boundary, q_boundary,
- project_to_boundary_first);
+ function, vec_result,
+ enforce_zero_boundary, q_boundary,
+ project_to_boundary_first);
}
template <int,int> class M_or_MC>
static inline
void do_interpolate_boundary_values (const M_or_MC<DH::dimension, DH::space_dimension> &,
- const DH &dof,
- const typename FunctionMap<DH::space_dimension>::type &function_map,
- std::map<types::global_dof_index,double> &boundary_values,
- const ComponentMask &component_mask,
- const dealii::internal::int2type<1>)
+ const DH &dof,
+ const typename FunctionMap<DH::space_dimension>::type &function_map,
+ std::map<types::global_dof_index,double> &boundary_values,
+ const ComponentMask &component_mask,
+ const dealii::internal::int2type<1>)
{
const unsigned int dim = DH::dimension;
const unsigned int spacedim=DH::space_dimension;
static inline
void
do_interpolate_boundary_values (const M_or_MC<DH::dimension, DH::space_dimension> &mapping,
- const DH &dof,
- const typename FunctionMap<DH::space_dimension>::type &function_map,
- std::map<types::global_dof_index,double> &boundary_values,
- const ComponentMask &component_mask,
- const dealii::internal::int2type<DH::dimension>)
+ const DH &dof,
+ const typename FunctionMap<DH::space_dimension>::type &function_map,
+ std::map<types::global_dof_index,double> &boundary_values,
+ const ComponentMask &component_mask,
+ const dealii::internal::int2type<DH::dimension>)
{
const unsigned int dim = DH::dimension;
const unsigned int spacedim=DH::space_dimension;
:
(dim == 2 ?
(i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+ :
+ (dim == 3 ?
+ (i<4*fe.dofs_per_vertex ?
+ i
:
- (dim == 3 ?
- (i<4*fe.dofs_per_vertex ?
- i
- :
- (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
- i+4*fe.dofs_per_vertex
- :
- i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+ (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
+ i+4*fe.dofs_per_vertex
:
- numbers::invalid_unsigned_int)));
+ i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+ :
+ numbers::invalid_unsigned_int)));
Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
// make sure that if this is not a primitive
}
else
// fe has only one component, so save some computations
- {
+ {
// get only the one component that this function has
dof_values_scalar.resize (fe.dofs_per_face);
function_map.find(boundary_component)->second
const ComponentMask &component_mask_)
{
do_interpolate_boundary_values (mapping, dof, function_map, boundary_values,
- component_mask_,
- dealii::internal::int2type<DH::dimension>());
+ component_mask_,
+ dealii::internal::int2type<DH::dimension>());
}
const ComponentMask &component_mask_)
{
do_interpolate_boundary_values (mapping, dof, function_map, boundary_values,
- component_mask_,
- dealii::internal::int2type<dim>());
+ component_mask_,
+ dealii::internal::int2type<dim>());
}
template <int> class Q_or_QC>
void
do_project_boundary_values (const M_or_MC<dim, spacedim> &mapping,
- const DH<dim, spacedim> &dof,
- const typename FunctionMap<spacedim>::type &boundary_functions,
- const Q_or_QC<dim-1> &q,
- std::map<types::global_dof_index,double> &boundary_values,
- std::vector<unsigned int> component_mapping)
+ const DH<dim, spacedim> &dof,
+ const typename FunctionMap<spacedim>::type &boundary_functions,
+ const Q_or_QC<dim-1> &q,
+ std::map<types::global_dof_index,double> &boundary_values,
+ std::vector<unsigned int> component_mapping)
{
// in 1d, projection onto the 0d end points == interpolation
if (dim == 1)
return;
}
- //TODO:[?] In project_boundary_values, no condensation of sparsity
- // structures, matrices and right hand sides or distribution of
- // solution vectors is performed. This is ok for dim<3 because then
- // there are no constrained nodes on the boundary, but is not
- // acceptable for higher dimensions. Fix this.
+ //TODO:[?] In project_boundary_values, no condensation of sparsity
+ // structures, matrices and right hand sides or distribution of
+ // solution vectors is performed. This is ok for dim<3 because then
+ // there are no constrained nodes on the boundary, but is not
+ // acceptable for higher dimensions. Fix this.
if (component_mapping.size() == 0)
{
// but it needs to be implemented
if (dim>=3)
{
- #ifdef DEBUG
- // Assert that there are no hanging nodes at the boundary
+#ifdef DEBUG
+ // Assert that there are no hanging nodes at the boundary
int level = -1;
for (typename DH<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
cell != dof.end(); ++cell)
}
}
}
- #endif
+#endif
}
sparsity.compress();
std::vector<unsigned int> component_mapping)
{
do_project_boundary_values (mapping, dof,
- boundary_functions,
- q, boundary_values,
- component_mapping);
+ boundary_functions,
+ q, boundary_values,
+ component_mapping);
}
std::vector<unsigned int> component_mapping)
{
project_boundary_values (hp::StaticMappingQ1<dim,spacedim>::mapping_collection, dof,
- boundary_function,
- q, boundary_values,
- component_mapping);
+ boundary_function,
+ q, boundary_values,
+ component_mapping);
}
template <int dim>
- std::ostream & operator << (std::ostream &out,
- const VectorDoFTuple<dim> &vdt)
+ std::ostream &operator << (std::ostream &out,
+ const VectorDoFTuple<dim> &vdt)
{
for (unsigned int d=0; d<dim; ++d)
- out << vdt.dof_indices[d] << (d < dim-1 ? " " : "");
+ out << vdt.dof_indices[d] << (d < dim-1 ? " " : "");
return out;
}
-= tol;
tangentials[q_point]
= (fe_values.get_mapping ()
- .transform_unit_to_real_cell (cell,
- shifted_reference_point_1)
- -
- fe_values.get_mapping ()
- .transform_unit_to_real_cell (cell,
- shifted_reference_point_2))
+ .transform_unit_to_real_cell (cell,
+ shifted_reference_point_1)
+ -
+ fe_values.get_mapping ()
+ .transform_unit_to_real_cell (cell,
+ shifted_reference_point_2))
/ tol;
tangentials[q_point]
/= std::sqrt (tangentials[q_point].square ());
fe_collection,
face_quadrature_collection,
update_q_points |
- update_normal_vectors);
+ update_normal_vectors);
// 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
vector_dofs.dof_indices[0] = face_dofs[i];
Assert(first_vector_component+dim<=fe.n_components(),
- ExcMessage("Error: the finite element does not have enough components "
- "to define a normal direction."));
+ ExcMessage("Error: the finite element does not have enough components "
+ "to define a normal direction."));
for (unsigned int k=0; k<fe.dofs_per_face; ++k)
if ((k != i)
// to check that they are small *relative to something
// else*). we do this and then normalize the length of the
// vector back to one, just to be on the safe side
- //
- // one more point: we would like to use the "real" normal
- // vector here, as provided by the boundary description
- // and as opposed to what we get from the FEValues object.
- // we do this in the immediately next line, but as is
- // obvious, the boundary only has a vague idea which side
- // of a cell it is on -- indicated by the face number. in
- // other words, it may provide the inner or outer normal.
- // by and large, there is no harm from this, since the
- // tangential vector we compute is still the same. however,
- // we do average over normal vectors from adjacent cells
- // and if they have recorded normal vectors from the inside
- // once and from the outside the other time, then this
- // averaging is going to run into trouble. as a consequence
- // we ask the mapping after all for its normal vector,
- // but we only ask it so that we can possibly correct the
- // sign of the normal vector provided by the boundary
- // if they should point in different directions. this is the
- // case in tests/deal.II/no_flux_11.
+ //
+ // one more point: we would like to use the "real" normal
+ // vector here, as provided by the boundary description
+ // and as opposed to what we get from the FEValues object.
+ // we do this in the immediately next line, but as is
+ // obvious, the boundary only has a vague idea which side
+ // of a cell it is on -- indicated by the face number. in
+ // other words, it may provide the inner or outer normal.
+ // by and large, there is no harm from this, since the
+ // tangential vector we compute is still the same. however,
+ // we do average over normal vectors from adjacent cells
+ // and if they have recorded normal vectors from the inside
+ // once and from the outside the other time, then this
+ // averaging is going to run into trouble. as a consequence
+ // we ask the mapping after all for its normal vector,
+ // but we only ask it so that we can possibly correct the
+ // sign of the normal vector provided by the boundary
+ // if they should point in different directions. this is the
+ // case in tests/deal.II/no_flux_11.
Point<dim> normal_vector
= (cell->face(face_no)->get_boundary()
.normal_vector (cell->face(face_no),
fe_values.quadrature_point(i)));
- if (normal_vector * fe_values.normal_vector(i) < 0)
- normal_vector *= -1;
+ if (normal_vector * fe_values.normal_vector(i) < 0)
+ normal_vector *= -1;
Assert (std::fabs(normal_vector.norm() - 1) < 1e-14,
ExcInternalError());
for (unsigned int d=0; d<dim; ++d)
std::make_pair (normal_vector,
cell)));
#ifdef DEBUG_NO_NORMAL_FLUX
- std::cout << "Adding normal vector:" << std::endl
- << " dofs=" << vector_dofs << std::endl
- << " cell=" << cell << " at " << cell->center() << std::endl
- << " normal=" << normal_vector << std::endl;
+ std::cout << "Adding normal vector:" << std::endl
+ << " dofs=" << vector_dofs << std::endl
+ << " cell=" << cell << " at " << cell->center() << std::endl
+ << " normal=" << normal_vector << std::endl;
#endif
}
}
same_dof_range[1] = dof_to_normals_map.end();
#ifdef DEBUG_NO_NORMAL_FLUX
- std::cout << "For dof indices <" << p->first << ">, found the following normals"
- << std::endl;
+ std::cout << "For dof indices <" << p->first << ">, found the following normals"
+ << std::endl;
for (typename DoFToNormalsMap::const_iterator
q = same_dof_range[0];
q != same_dof_range[1]; ++q)
- std::cout << " " << q->second.first
- << " from cell " << q->second.second
- << std::endl;
+ std::cout << " " << q->second.first
+ << " from cell " << q->second.second
+ << std::endl;
#endif
Assert (cell_to_normals_map.size() >= 1, ExcInternalError());
#ifdef DEBUG_NO_NORMAL_FLUX
- std::cout << " cell_to_normals_map:" << std::endl;
+ std::cout << " cell_to_normals_map:" << std::endl;
for (typename CellToNormalsMap::const_iterator
x = cell_to_normals_map.begin();
x != cell_to_normals_map.end(); ++x)
- std::cout << " " << x->first << " -> ("
- << x->second.first << ',' << x->second.second << ')'
- << std::endl;
+ std::cout << " " << x->first << " -> ("
+ << x->second.first << ',' << x->second.second << ')'
+ << std::endl;
#endif
// count the maximum number of contributions from each cell
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
- //
- // the typical case is in 2d where multiple cells meet at one
- // vertex sitting on the boundary. same in 3d for a vertex that
- // is associated with only one of the boundary indicators passed
- // to this function
+ // 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
+ //
+ // the typical case is in 2d where multiple cells meet at one
+ // vertex sitting on the boundary. same in 3d for a vertex that
+ // is associated with only one of the boundary indicators passed
+ // to this function
case 1:
{
// if this is dim-1 for all cells (if a single cell has
// contributed dim, or if all adjacent cells have contributed 1
// normal vector, this is already handled above).
- //
- // we only implement the case that all cells contribute
- // dim-1 because we assume that we are following an edge
- // of the domain (think: we are looking at a vertex
- // located on one of the edges of a refined cube where the
- // boundary indicators of the two adjacent faces of the
- // cube are both listed in the set of boundary indicators
- // passed to this function). in that case, all cells along
- // that edge of the domain are assumed to have contributed
- // dim-1 normal vectors. however, there are cases where
- // this assumption is not justified (see the lengthy
- // explanation in test no_flux_12.cc) and in those cases
- // we simply ignore the cell that contributes only
- // once. this is also discussed at length in the
- // documentation of this function.
+ //
+ // we only implement the case that all cells contribute
+ // dim-1 because we assume that we are following an edge
+ // of the domain (think: we are looking at a vertex
+ // located on one of the edges of a refined cube where the
+ // boundary indicators of the two adjacent faces of the
+ // cube are both listed in the set of boundary indicators
+ // passed to this function). in that case, all cells along
+ // that edge of the domain are assumed to have contributed
+ // dim-1 normal vectors. however, there are cases where
+ // this assumption is not justified (see the lengthy
+ // explanation in test no_flux_12.cc) and in those cases
+ // we simply ignore the cell that contributes only
+ // once. this is also discussed at length in the
+ // documentation of this function.
//
// for each contributing cell compute the tangential vector that
// remains unconstrained
++contribution)
{
#ifdef DEBUG_NO_NORMAL_FLUX
- std::cout << " Treating edge case with dim-1 contributions." << std::endl
- << " Looking at cell " << contribution->first
- << " which has contributed these normal vectors:"
- << std::endl;
- for (typename std::list<Tensor<1,dim> >::const_iterator
- t = contribution->second.begin();
- t != contribution->second.end();
- ++t)
- std::cout << " " << *t << std::endl;
+ std::cout << " Treating edge case with dim-1 contributions." << std::endl
+ << " Looking at cell " << contribution->first
+ << " which has contributed these normal vectors:"
+ << std::endl;
+ for (typename std::list<Tensor<1,dim> >::const_iterator
+ t = contribution->second.begin();
+ t != contribution->second.end();
+ ++t)
+ std::cout << " " << *t << std::endl;
#endif
- // as mentioned above, simply ignore cells that only
- // contribute once
+ // as mentioned above, simply ignore cells that only
+ // contribute once
if (contribution->second.size() < dim-1)
- continue;
+ continue;
Tensor<1,dim> normals[dim-1];
{
template <int dim, template <int, int> class DH, int spacedim>
void
- compute_normal_flux_constraints (const DH<dim,spacedim>&dof_handler,
+ compute_normal_flux_constraints (const DH<dim,spacedim> &dof_handler,
const unsigned int first_vector_component,
const std::set<types::boundary_id> &boundary_ids,
ConstraintMatrix &constraints,
std::max(n_scalar_indices,
fe.face_system_to_component_index(i).second+1);
cell_vector_dofs[fe.face_system_to_component_index(i).second]
- [fe.face_system_to_component_index(i).first-first_vector_component]
+ [fe.face_system_to_component_index(i).first-first_vector_component]
= face_dofs[i];
}
// can find constrained ones
unsigned int n_total_constraints_found = 0;
for (typename std::set<std_cxx1x::array<types::global_dof_index,dim>,
- PointComparator<dim> >::const_iterator it=vector_dofs.begin();
- it!=vector_dofs.end(); ++it)
+ PointComparator<dim> >::const_iterator it=vector_dofs.begin();
+ it!=vector_dofs.end(); ++it)
{
unsigned int n_constraints = 0;
bool is_constrained[dim];
normal[d] = 1.;
}
AssertIndexRange(constrained_index, dim);
- const std::vector<std::pair<types::global_dof_index, double> >* constrained
+ const std::vector<std::pair<types::global_dof_index, double> > *constrained
= no_normal_flux_constraints.get_constraint_entries((*it)[constrained_index]);
// find components to which this index is constrained to
Assert(constrained != 0, ExcInternalError());
double sum = 0;
for (unsigned int k=0; k<n_components; ++k)
sum += data.psi_values[q](k) * data.psi_values[q](k) *
- data.weight_vectors[q](k);
+ data.weight_vectors[q](k);
diff += sum * fe_values.JxW(q);
}
// Compute the root only, if no derivative values are added later
double sum = 0;
for (unsigned int k=0; k<n_components; ++k)
sum += (data.psi_grads[q][k] * data.psi_grads[q][k]) *
- data.weight_vectors[q](k);
+ data.weight_vectors[q](k);
diff += sum * fe_values.JxW(q);
}
diff = std::sqrt(diff);
case W1infty_seminorm:
case W1infty_norm:
- {
- double t = 0;
- for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int k=0; k<n_components; ++k)
- for (unsigned int d=0; d<dim; ++d)
- t = std::max(t, std::abs(data.psi_grads[q][k][d]) *
- data.weight_vectors[q](k));
+ {
+ double t = 0;
+ for (unsigned int q=0; q<n_q_points; ++q)
+ for (unsigned int k=0; k<n_components; ++k)
+ for (unsigned int d=0; d<dim; ++d)
+ t = std::max(t, std::abs(data.psi_grads[q][k][d]) *
+ data.weight_vectors[q](k));
- // then add seminorm to norm if that had previously been computed
- diff += t;
- }
- break;
+ // then add seminorm to norm if that had previously been computed
+ diff += t;
+ }
+ break;
default:
break;
}
case W1infty_seminorm:
update_flags |= UpdateFlags (update_gradients);
if (spacedim == dim+1)
- update_flags |= UpdateFlags (update_normal_vectors);
+ update_flags |= UpdateFlags (update_normal_vectors);
break;
case W1infty_norm:
update_flags |= UpdateFlags (update_gradients);
if (spacedim == dim+1)
- update_flags |= UpdateFlags (update_normal_vectors);
- // no break!
+ update_flags |= UpdateFlags (update_normal_vectors);
+ // no break!
default:
update_flags |= UpdateFlags (update_values);
Point<dim> q1, q2, q3, q4;
Vector<double>
v1(this->n_components), v2(this->n_components),
- v3(this->n_components), v4(this->n_components);
+ v3(this->n_components), v4(this->n_components);
const double h_inv_12=1./(12*h);
for (unsigned int i=0; i<dim; ++i)
{
template <>
void DataOutBase::write_svg<1,1> (const std::vector<Patch<1,1> > &patches,
- const std::vector<std::string> &data_names,
- const std::vector<std_cxx1x::tuple<unsigned int, unsigned int, std::string> > &vector_data_ranges,
- const SvgFlags &flags,
- std::ostream &out)
+ const std::vector<std::string> &data_names,
+ const std::vector<std_cxx1x::tuple<unsigned int, unsigned int, std::string> > &vector_data_ranges,
+ const SvgFlags &flags,
+ std::ostream &out)
{
AssertThrow (false, ExcNotImplemented());
}
MPI_Info_create(&info);
MPI_File fh;
err = MPI_File_open(comm, const_cast<char *>(filename),
- MPI_MODE_CREATE | MPI_MODE_WRONLY, info, &fh);
+ MPI_MODE_CREATE | MPI_MODE_WRONLY, info, &fh);
AssertThrow(err==0, ExcMessage("Unable to open file with MPI_File_open!"));
#endif
}
-const char* ExceptionBase::what() const throw()
+const char *ExceptionBase::what() const throw()
{
// If no error c_string was generated so far, do it now:
if (what_str == "")
namespace fparser
{
class FunctionParser: public ::FunctionParser
-{};
+ {};
}
#else
namespace fparser
{
class FunctionParser
-{};
+ {};
}
#endif
GeometryInfo<dim>::
alternating_form_at_vertices
#ifndef DEAL_II_CONSTEXPR_BUG
- (const Point<spacedim> (&vertices)[vertices_per_cell],
- Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell])
+(const Point<spacedim> (&vertices)[vertices_per_cell],
+ Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell])
#else
- (const Point<spacedim> *vertices,
- Tensor<spacedim-dim,spacedim> *forms)
+(const Point<spacedim> *vertices,
+ Tensor<spacedim-dim,spacedim> *forms)
#endif
{
// for each of the vertices,
{
// Implement a minimalistic stream buffer that only stores the fact
// whether overflow or sync was called
- public:
- QueryStreambuf()
- : flushed_(false), newline_written_(false)
- {
- }
- bool flushed() { return flushed_; }
- bool newline_written() { return newline_written_; }
- private:
- int_type overflow(int_type ch)
- {
- newline_written_ = true;
- return ch;
- }
- int sync()
- {
- flushed_ = true;
- return 0;
- }
- bool flushed_;
- bool newline_written_;
+ public:
+ QueryStreambuf()
+ : flushed_(false), newline_written_(false)
+ {
+ }
+ bool flushed()
+ {
+ return flushed_;
+ }
+ bool newline_written()
+ {
+ return newline_written_;
+ }
+ private:
+ int_type overflow(int_type ch)
+ {
+ newline_written_ = true;
+ return ch;
+ }
+ int sync()
+ {
+ flushed_ = true;
+ return 0;
+ }
+ bool flushed_;
+ bool newline_written_;
} query_streambuf;
{
void max_reduce ( const void *in_lhs_,
void *inout_rhs_,
int *len,
- MPI_Datatype * )
+ MPI_Datatype *)
{
const MinMaxAvg *in_lhs = static_cast<const MinMaxAvg *>(in_lhs_);
MinMaxAvg *inout_rhs = static_cast<MinMaxAvg *>(inout_rhs_);
// when running PETSc, because we initialize MPI ourselves before calling
// PetscInitialize
#ifdef DEAL_II_WITH_MPI
- int mpi_err = 0;
+ int mpi_err = 0;
int MPI_has_been_started = 0;
MPI_Initialized(&MPI_has_been_started);
void MultithreadInfo::set_thread_limit(const unsigned int max_threads)
{
Assert(n_max_threads==numbers::invalid_unsigned_int,
- ExcMessage("Calling set_thread_limit() more than once is not supported!"));
+ ExcMessage("Calling set_thread_limit() more than once is not supported!"));
unsigned int max_threads_env = numbers::invalid_unsigned_int;
char *penv;
// the whole filename specified, including (possibly)
// the suffix
{
- real_name = *path + filename;
- if (debug > 1)
- deallog << "PathSearch[" << cls << "] trying "
- << real_name << std::endl;
- FILE *fp = fopen(real_name.c_str(), open_mode);
- if (fp != 0)
- {
- if (debug > 0)
- deallog << "PathSearch[" << cls << "] opened "
- << real_name << std::endl;
- fclose(fp);
- return real_name;
- }
+ real_name = *path + filename;
+ if (debug > 1)
+ deallog << "PathSearch[" << cls << "] trying "
+ << real_name << std::endl;
+ FILE *fp = fopen(real_name.c_str(), open_mode);
+ if (fp != 0)
+ {
+ if (debug > 0)
+ deallog << "PathSearch[" << cls << "] opened "
+ << real_name << std::endl;
+ fclose(fp);
+ return real_name;
+ }
}
// try again with the suffix appended, unless there is
real_name = *path + filename + suffix;
if (debug > 1)
deallog << "PathSearch[" << cls << "] trying "
- << real_name << std::endl;
+ << real_name << std::endl;
FILE *fp = fopen(real_name.c_str(), open_mode);
if (fp != 0)
{
if (debug > 0)
deallog << "PathSearch[" << cls << "] opened "
- << real_name << std::endl;
+ << real_name << std::endl;
fclose(fp);
return real_name;
}
}
break;
- // manually implement size 1 (values only), size 2 (value + first
- // derivative), and size 3 (up to second derivative) since they
- // might be called often. then, we can unroll the loop.
+ // manually implement size 1 (values only), size 2 (value + first
+ // derivative), and size 3 (up to second derivative) since they
+ // might be called often. then, we can unroll the loop.
case 1:
values[0] = 1;
for (unsigned int i=0; i<n_supp; ++i)
template <int dim>
Quadrature<dim>::Quadrature (const Quadrature<dim != 1 ? 1 : 0> &q)
- :
- Subscriptor(),
- quadrature_points (Utilities::fixed_power<dim>(q.size())),
- weights (Utilities::fixed_power<dim>(q.size()))
+:
+Subscriptor(),
+ quadrature_points (Utilities::fixed_power<dim>(q.size())),
+ weights (Utilities::fixed_power<dim>(q.size()))
{
Assert (dim <= 3, ExcNotImplemented());
out << std::flush;
return;
}
-
+
case simple_table_with_separate_column_description:
{
// write the captions
out << std::flush;
return;
}
-
+
case table_with_separate_column_description:
{
// writing the captions for table_with_separate_column_description
}
break;
}
-
+
case table_with_headers:
{
// This format output supercolumn headers and aligns them centered
out << '\n';
break;
}
-
+
default:
Assert (false, ExcInternalError());
}
-
+
// finally output the data itself for
// table_with_headers or table_with_separate_column_description:
for (unsigned int i=0; i<nrows; ++i)
int face_left,
int face_right,
int orientation)
- = p4est_connectivity_join_faces;
+ = p4est_connectivity_join_faces;
#endif
void (&functions<2>::connectivity_destroy) (p4est_connectivity_t *connectivity)
= p4est_connectivity_save;
int (&functions<2>::connectivity_is_valid) (types<2>::connectivity
- *connectivity)
+ *connectivity)
= p4est_connectivity_is_valid;
types<2>::connectivity *
= p8est_connectivity_save;
int (&functions<3>::connectivity_is_valid) (types<3>::connectivity
- *connectivity)
+ *connectivity)
= p8est_connectivity_is_valid;
types<3>::connectivity *
std::vector<unsigned int> &vertex_touch_count,
std::vector<std::list<
std::pair<typename Triangulation<dim,spacedim>::active_cell_iterator,unsigned int> > >
- & vertex_to_cell)
+ &vertex_to_cell)
{
vertex_touch_count.resize (triangulation.n_vertices());
vertex_to_cell.resize (triangulation.n_vertices());
std::vector<unsigned int> &edge_touch_count,
std::vector<std::list<
std::pair<typename Triangulation<dim,spacedim>::active_cell_iterator,unsigned int> > >
- & edge_to_cell)
+ &edge_to_cell)
{
Assert (triangulation.n_levels() == 1, ExcInternalError());
const std::vector<unsigned int> &vertex_touch_count,
const std::vector<std::list<
std::pair<typename Triangulation<dim,spacedim>::active_cell_iterator,unsigned int> > >
- & vertex_to_cell,
+ &vertex_to_cell,
const std::vector<types::global_dof_index> &coarse_cell_to_p4est_tree_permutation,
const bool set_vertex_info,
typename internal::p4est::types<dim>::connectivity *connectivity)
Triangulation<dim,spacedim>::
load(const char *filename)
{
- if (parallel_ghost != 0) {
- dealii::internal::p4est::functions<dim>::ghost_destroy (parallel_ghost);
- parallel_ghost = 0;
- }
+ if (parallel_ghost != 0)
+ {
+ dealii::internal::p4est::functions<dim>::ghost_destroy (parallel_ghost);
+ parallel_ghost = 0;
+ }
dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
parallel_forest = 0;
dealii::internal::p4est::functions<dim>::connectivity_destroy (connectivity);
// this command changes the distribution for some reason, so we
// will leave it in here.
dealii::internal::p4est::functions<dim>::
- partition (parallel_forest,
- /* prepare coarsening */ 1,
- /* weight_callback */ NULL);
+ partition (parallel_forest,
+ /* prepare coarsening */ 1,
+ /* weight_callback */ NULL);
}
num_ett = std::accumulate (edge_touch_count.begin(),
edge_touch_count.end(),
0);
-
+
// now create a connectivity object with the right sizes for all arrays
const bool set_vertex_info
#ifdef DEBUG
// function. Note that the smoothing flag is used in the normal
// refinement process.
typename Triangulation<dim,spacedim>::MeshSmoothing
- save_smooth = this->smooth_grid;
+ save_smooth = this->smooth_grid;
// We will refine manually to match the p4est further down, which
// obeys a level difference of 2 at each vertex (see the balance call
parallel_ghost = 0;
}
parallel_ghost = dealii::internal::p4est::functions<dim>::ghost_new (parallel_forest,
- (dim == 2
- ?
- typename dealii::internal::p4est::types<dim>::
- balance_type(P4EST_CONNECT_CORNER)
- :
- typename dealii::internal::p4est::types<dim>::
- balance_type(P8EST_CONNECT_CORNER)));
+ (dim == 2
+ ?
+ typename dealii::internal::p4est::types<dim>::
+ balance_type(P4EST_CONNECT_CORNER)
+ :
+ typename dealii::internal::p4est::types<dim>::
+ balance_type(P8EST_CONNECT_CORNER)));
Assert (parallel_ghost, ExcInternalError());
Assert(cell->neighbor(f)->active(), ExcInternalError());
Assert(cell->neighbor(f)->subdomain_id() != numbers::artificial_subdomain_id, ExcInternalError());
Assert(cell->neighbor(f)->level_subdomain_id() == numbers::artificial_subdomain_id
- || cell->neighbor(f)->level_subdomain_id() == cell->neighbor(f)->subdomain_id(), ExcInternalError());
+ || cell->neighbor(f)->level_subdomain_id() == cell->neighbor(f)->subdomain_id(), ExcInternalError());
cell->neighbor(f)->set_level_subdomain_id(cell->neighbor(f)->subdomain_id());
}
}
register_data_attach (const std::size_t size,
const std_cxx1x::function<void(const cell_iterator &,
const CellStatus,
- void *)> & pack_callback)
+ void *)> &pack_callback)
{
Assert(size>0, ExcMessage("register_data_attach(), size==0"));
Assert(attached_data_pack_callbacks.size()==n_attached_datas,
notify_ready_to_unpack (const unsigned int offset,
const std_cxx1x::function<void (const cell_iterator &,
const CellStatus,
- const void *)> & unpack_callback)
+ const void *)> &unpack_callback)
{
Assert (offset < attached_data_size, ExcMessage ("invalid offset in notify_ready_to_unpack()"));
Assert (n_attached_datas > 0, ExcMessage ("notify_ready_to_unpack() called too often"));
struct find_ghosts
{
typename dealii::parallel::distributed::Triangulation<dim,spacedim> *triangulation;
- sc_array_t * subids;
+ sc_array_t *subids;
std::map<unsigned int, std::set<dealii::types::subdomain_id> >
*vertices_with_ghost_neighbors;
};
template <int dim, int spacedim>
void
find_ghosts_corner
- (typename dealii::internal::p4est::iter<dim>::corner_info * info,
+ (typename dealii::internal::p4est::iter<dim>::corner_info *info,
void *user_data)
{
int i, j;
int nsides = info->sides.elem_count;
- typename dealii::internal::p4est::iter<dim>::corner_side * sides =
+ typename dealii::internal::p4est::iter<dim>::corner_side *sides =
(typename dealii::internal::p4est::iter<dim>::corner_side *)
(info->sides.array);
struct find_ghosts<dim,spacedim> *fg = static_cast<struct find_ghosts<dim,spacedim> *>(user_data);
*vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors;
subids->elem_count = 0;
- for (i = 0; i < nsides; i++) {
- if (sides[i].is_ghost) {
- typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad));
- Assert (cell->is_ghost(), ExcMessage ("ghost quad did not find ghost cell"));
- dealii::types::subdomain_id *subid =
- static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
- *subid = cell->subdomain_id();
+ for (i = 0; i < nsides; i++)
+ {
+ if (sides[i].is_ghost)
+ {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad));
+ Assert (cell->is_ghost(), ExcMessage ("ghost quad did not find ghost cell"));
+ dealii::types::subdomain_id *subid =
+ static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
+ *subid = cell->subdomain_id();
+ }
}
- }
- if (!subids->elem_count) {
- return;
- }
+ if (!subids->elem_count)
+ {
+ return;
+ }
nsubs = (int) subids->elem_count;
subdomain_ids = (dealii::types::subdomain_id *) (subids->array);
- for (i = 0; i < nsides; i++) {
- if (!sides[i].is_ghost) {
- typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad));
+ for (i = 0; i < nsides; i++)
+ {
+ if (!sides[i].is_ghost)
+ {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].quad));
- Assert (!cell->is_ghost(), ExcMessage ("local quad found ghost cell"));
+ Assert (!cell->is_ghost(), ExcMessage ("local quad found ghost cell"));
- for (j = 0; j < nsubs; j++) {
- (*vertices_with_ghost_neighbors)[cell->vertex_index(sides[i].corner)]
- .insert (subdomain_ids[j]);
- }
+ for (j = 0; j < nsubs; j++)
+ {
+ (*vertices_with_ghost_neighbors)[cell->vertex_index(sides[i].corner)]
+ .insert (subdomain_ids[j]);
+ }
+ }
}
- }
subids->elem_count = 0;
}
template <int dim, int spacedim>
void
find_ghosts_edge
- (typename dealii::internal::p4est::iter<dim>::edge_info * info,
+ (typename dealii::internal::p4est::iter<dim>::edge_info *info,
void *user_data)
{
int i, j, k;
int nsides = info->sides.elem_count;
- typename dealii::internal::p4est::iter<dim>::edge_side * sides =
+ typename dealii::internal::p4est::iter<dim>::edge_side *sides =
(typename dealii::internal::p4est::iter<dim>::edge_side *)
(info->sides.array);
struct find_ghosts<dim,spacedim> *fg = static_cast<struct find_ghosts<dim,spacedim> *>(user_data);
*vertices_with_ghost_neighbors = fg->vertices_with_ghost_neighbors;
subids->elem_count = 0;
- for (i = 0; i < nsides; i++) {
- if (sides[i].is_hanging) {
- for (j = 0; j < 2; j++) {
- if (sides[i].is.hanging.is_ghost[j]) {
- typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
- dealii::types::subdomain_id *subid =
- static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
- *subid = cell->subdomain_id();
+ for (i = 0; i < nsides; i++)
+ {
+ if (sides[i].is_hanging)
+ {
+ for (j = 0; j < 2; j++)
+ {
+ if (sides[i].is.hanging.is_ghost[j])
+ {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
+ dealii::types::subdomain_id *subid =
+ static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
+ *subid = cell->subdomain_id();
+ }
+ }
}
- }
}
- }
- if (!subids->elem_count) {
- return;
- }
+ if (!subids->elem_count)
+ {
+ return;
+ }
nsubs = (int) subids->elem_count;
subdomain_ids = (dealii::types::subdomain_id *) (subids->array);
- for (i = 0; i < nsides; i++) {
- if (sides[i].is_hanging) {
- for (j = 0; j < 2; j++) {
- if (!sides[i].is.hanging.is_ghost[j]) {
- typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
+ for (i = 0; i < nsides; i++)
+ {
+ if (sides[i].is_hanging)
+ {
+ for (j = 0; j < 2; j++)
+ {
+ if (!sides[i].is.hanging.is_ghost[j])
+ {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
- for (k = 0; k < nsubs; k++) {
- (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_edge_corners[sides[i].edge][1^j])]
- .insert (subdomain_ids[k]);
- }
+ for (k = 0; k < nsubs; k++)
+ {
+ (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_edge_corners[sides[i].edge][1^j])]
+ .insert (subdomain_ids[k]);
+ }
+ }
+ }
}
- }
}
- }
subids->elem_count = 0;
}
template <int dim, int spacedim>
void
find_ghosts_face
- (typename dealii::internal::p4est::iter<dim>::face_info * info,
+ (typename dealii::internal::p4est::iter<dim>::face_info *info,
void *user_data)
{
int i, j, k;
int nsides = info->sides.elem_count;
- typename dealii::internal::p4est::iter<dim>::face_side * sides =
+ typename dealii::internal::p4est::iter<dim>::face_side *sides =
(typename dealii::internal::p4est::iter<dim>::face_side *)
(info->sides.array);
struct find_ghosts<dim,spacedim> *fg = static_cast<struct find_ghosts<dim,spacedim> *>(user_data);
int limit = (dim == 2) ? 2 : 4;
subids->elem_count = 0;
- for (i = 0; i < nsides; i++) {
- if (sides[i].is_hanging) {
- for (j = 0; j < limit; j++) {
- if (sides[i].is.hanging.is_ghost[j]) {
- typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
- dealii::types::subdomain_id *subid =
- static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
- *subid = cell->subdomain_id();
+ for (i = 0; i < nsides; i++)
+ {
+ if (sides[i].is_hanging)
+ {
+ for (j = 0; j < limit; j++)
+ {
+ if (sides[i].is.hanging.is_ghost[j])
+ {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
+ dealii::types::subdomain_id *subid =
+ static_cast<dealii::types::subdomain_id *>(sc_array_push (subids));
+ *subid = cell->subdomain_id();
+ }
+ }
}
- }
}
- }
- if (!subids->elem_count) {
- return;
- }
+ if (!subids->elem_count)
+ {
+ return;
+ }
nsubs = (int) subids->elem_count;
subdomain_ids = (dealii::types::subdomain_id *) (subids->array);
- for (i = 0; i < nsides; i++) {
- if (sides[i].is_hanging) {
- for (j = 0; j < limit; j++) {
- if (!sides[i].is.hanging.is_ghost[j]) {
- typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
+ for (i = 0; i < nsides; i++)
+ {
+ if (sides[i].is_hanging)
+ {
+ for (j = 0; j < limit; j++)
+ {
+ if (!sides[i].is.hanging.is_ghost[j])
+ {
+ typename dealii::parallel::distributed::Triangulation<dim,spacedim>::cell_iterator cell = cell_from_quad (triangulation, sides[i].treeid, *(sides[i].is.hanging.quad[j]));
- for (k = 0; k < nsubs; k++) {
- if (dim == 2) {
- (*vertices_with_ghost_neighbors)[cell->vertex_index(p4est_face_corners[sides[i].face][(limit - 1)^j])]
- .insert (subdomain_ids[k]);
- }
- else {
- (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_face_corners[sides[i].face][(limit - 1)^j])]
- .insert (subdomain_ids[k]);
+ for (k = 0; k < nsubs; k++)
+ {
+ if (dim == 2)
+ {
+ (*vertices_with_ghost_neighbors)[cell->vertex_index(p4est_face_corners[sides[i].face][(limit - 1)^j])]
+ .insert (subdomain_ids[k]);
+ }
+ else
+ {
+ (*vertices_with_ghost_neighbors)[cell->vertex_index(p8est_face_corners[sides[i].face][(limit - 1)^j])]
+ .insert (subdomain_ids[k]);
+ }
+ }
+ }
}
- }
}
- }
}
- }
subids->elem_count = 0;
}
void
Triangulation<1,1>::
fill_vertices_with_ghost_neighbors
- (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- &vertices_with_ghost_neighbors)
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
{
Assert (false, ExcNotImplemented());
}
void
Triangulation<1,2>::
fill_vertices_with_ghost_neighbors
- (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- &vertices_with_ghost_neighbors)
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
{
Assert (false, ExcNotImplemented());
}
void
Triangulation<1,3>::
fill_vertices_with_ghost_neighbors
- (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- &vertices_with_ghost_neighbors)
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
{
Assert (false, ExcNotImplemented());
}
void
Triangulation<2,2>::
fill_vertices_with_ghost_neighbors
- (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- &vertices_with_ghost_neighbors)
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
{
struct find_ghosts<2,2> fg;
void
Triangulation<2,3>::
fill_vertices_with_ghost_neighbors
- (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- &vertices_with_ghost_neighbors)
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
{
struct find_ghosts<2,3> fg;
void
Triangulation<3,3>::
fill_vertices_with_ghost_neighbors
- (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
- &vertices_with_ghost_neighbors)
+ (std::map<unsigned int, std::set<dealii::types::subdomain_id> >
+ &vertices_with_ghost_neighbors)
{
struct find_ghosts<3,3> fg;
template<int dim, int spacedim>
void
Triangulation<dim,spacedim>::add_periodicity
- (const std::vector<GridTools::PeriodicFacePair<cell_iterator> >&
- periodicity_vector)
+ (const std::vector<GridTools::PeriodicFacePair<cell_iterator> > &
+ periodicity_vector)
{
- #if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
+#if DEAL_II_P4EST_VERSION_GTE(0,3,4,1)
Assert (triangulation_has_content == true,
ExcMessage ("The triangulation is empty!"));
Assert (this->n_levels() == 1,
ExcMessage ("The triangulation is refined!"));
-
+
typedef std::vector<GridTools::PeriodicFacePair<cell_iterator> >
- FaceVector;
+ FaceVector;
typename FaceVector::const_iterator it, periodic_end;
it = periodicity_vector.begin();
periodic_end = periodicity_vector.end();
-
+
for (; it<periodic_end; ++it)
- {
- const cell_iterator first_cell = it->cell[0];
- const cell_iterator second_cell = it->cell[1];
- const unsigned int face_left = it->face_idx[0];
- const unsigned int face_right = it->face_idx[1];
-
- //respective cells of the matching faces in p4est
- const unsigned int tree_left
- = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
- first_cell)];
- const unsigned int tree_right
- = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
- second_cell)];
-
- //TODO Add support for non default orientation.
- Assert(it->orientation == 1,
- ExcMessage("Found a face match with non standard orientation. "
- "This function is only suitable for meshes with "
- "cells in default orientation"));
-
- dealii::internal::p4est::functions<dim>::
+ {
+ const cell_iterator first_cell = it->cell[0];
+ const cell_iterator second_cell = it->cell[1];
+ const unsigned int face_left = it->face_idx[0];
+ const unsigned int face_right = it->face_idx[1];
+
+ //respective cells of the matching faces in p4est
+ const unsigned int tree_left
+ = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
+ first_cell)];
+ const unsigned int tree_right
+ = coarse_cell_to_p4est_tree_permutation[std::distance(this->begin(),
+ second_cell)];
+
+ //TODO Add support for non default orientation.
+ Assert(it->orientation == 1,
+ ExcMessage("Found a face match with non standard orientation. "
+ "This function is only suitable for meshes with "
+ "cells in default orientation"));
+
+ dealii::internal::p4est::functions<dim>::
connectivity_join_faces (connectivity,
tree_left,
tree_right,
face_right,
/* orientation */ 0);
- /* The orientation parameter above describes the difference between
- * the cell on the left and the cell on the right would number of the
- * corners of the face. In the periodic domains most users will want,
- * this orientation will be 0, i.e. the two cells would number the
- * corners the same way. More exotic periodicity, like moebius strips
- * or converting an unstructured quad/hex mesh into a periodic domain,
- * are not supported right now, and undefined behavior will occur if
- * users try to make them periodic. This may be addressed at a later
- * date.
- */
- }
-
-
+ /* The orientation parameter above describes the difference between
+ * the cell on the left and the cell on the right would number of the
+ * corners of the face. In the periodic domains most users will want,
+ * this orientation will be 0, i.e. the two cells would number the
+ * corners the same way. More exotic periodicity, like moebius strips
+ * or converting an unstructured quad/hex mesh into a periodic domain,
+ * are not supported right now, and undefined behavior will occur if
+ * users try to make them periodic. This may be addressed at a later
+ * date.
+ */
+ }
+
+
Assert(dealii::internal::p4est::functions<dim>::connectivity_is_valid
- (connectivity) == 1, ExcInternalError());
-
+ (connectivity) == 1, ExcInternalError());
+
// now create a forest out of the connectivity data structure
dealii::internal::p4est::functions<dim>::destroy (parallel_forest);
parallel_forest
= dealii::internal::p4est::functions<dim>::
- new_forest (mpi_communicator,
- connectivity,
- /* minimum initial number of quadrants per tree */ 0,
- /* minimum level of upfront refinement */ 0,
- /* use uniform upfront refinement */ 1,
- /* user_data_size = */ 0,
- /* user_data_constructor = */ NULL,
- /* user_pointer */ this);
+ new_forest (mpi_communicator,
+ connectivity,
+ /* minimum initial number of quadrants per tree */ 0,
+ /* minimum level of upfront refinement */ 0,
+ /* use uniform upfront refinement */ 1,
+ /* user_data_size = */ 0,
+ /* user_data_constructor = */ NULL,
+ /* user_pointer */ this);
try
- {
- copy_local_forest_to_triangulation ();
- }
+ {
+ copy_local_forest_to_triangulation ();
+ }
catch (const typename Triangulation<dim>::DistortedCellList &)
- {
- // the underlying triangulation should not be checking for distorted
- // cells
- AssertThrow (false, ExcInternalError());
- }
-
- #else
+ {
+ // the underlying triangulation should not be checking for distorted
+ // cells
+ AssertThrow (false, ExcInternalError());
+ }
+
+#else
Assert(false, ExcMessage ("Need p4est version >= 0.3.4.1!"));
- #endif
+#endif
}
const std::vector<types::global_dof_index> perm = triangulation->get_p4est_tree_to_coarse_cell_permutation ();
types::global_dof_index dealii_index = perm[treeidx];
- for (i = 0; i < l; i++) {
- typename dealii::Triangulation<dim,spacedim>::cell_iterator cell (triangulation, i, dealii_index);
- child_id = internal::p4est::functions<dim>::quadrant_ancestor_id (&quad, i + 1);
- Assert (cell->has_children (), ExcMessage ("p4est quadrant does not correspond to a cell!"));
- dealii_index = cell->child_index(child_id);
- }
+ for (i = 0; i < l; i++)
+ {
+ typename dealii::Triangulation<dim,spacedim>::cell_iterator cell (triangulation, i, dealii_index);
+ child_id = internal::p4est::functions<dim>::quadrant_ancestor_id (&quad, i + 1);
+ Assert (cell->has_children (), ExcMessage ("p4est quadrant does not correspond to a cell!"));
+ dealii_index = cell->child_index(child_id);
+ }
typename dealii::Triangulation<dim,spacedim>::cell_iterator out_cell (triangulation, l, dealii_index);
namespace internal
{
template<int dim, int spacedim>
- std::string policy_to_string(const dealii::internal::DoFHandler::Policy::PolicyBase<dim,spacedim> & policy)
+ std::string policy_to_string(const dealii::internal::DoFHandler::Policy::PolicyBase<dim,spacedim> &policy)
{
std::string policy_name;
if (dynamic_cast<const typename dealii::internal::DoFHandler::Policy::Sequential<dim,spacedim>*>(&policy))
- policy_name = "Policy::Sequential<";
+ policy_name = "Policy::Sequential<";
+ else if (dynamic_cast<const typename dealii::internal::DoFHandler::Policy::ParallelDistributed<dim,spacedim>*>(&policy))
+ policy_name = "Policy::ParallelDistributed<";
else
- if (dynamic_cast<const typename dealii::internal::DoFHandler::Policy::ParallelDistributed<dim,spacedim>*>(&policy))
- policy_name = "Policy::ParallelDistributed<";
- else
- AssertThrow(false, ExcNotImplemented());
+ AssertThrow(false, ExcNotImplemented());
policy_name += Utilities::int_to_string(dim)+
- ","+Utilities::int_to_string(spacedim)+">";
+ ","+Utilities::int_to_string(spacedim)+">";
return policy_name;
}
16*dof_handler.selected_fe->dofs_per_quad;
break;
- // the following numbers are not based on actual counting but by
- // extrapolating the number sequences from the previous ones (for
- // example, for dofs_per_vertex, the sequence above is 19, 21, 28,
- // 30, 37, and is continued as follows):
+ // the following numbers are not based on actual counting but by
+ // extrapolating the number sequences from the previous ones (for
+ // example, for dofs_per_vertex, the sequence above is 19, 21, 28,
+ // 30, 37, and is continued as follows):
case 9:
max_couplings=39*dof_handler.selected_fe->dofs_per_vertex +
59*dof_handler.selected_fe->dofs_per_line +
//have us as a neighbor. In this case we and the children are active.
if (dealii_cell->active())
{
- for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
if (dealii_cell->at_boundary(f))
continue;
:
(dim == 2 ?
(i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+ :
+ (dim == 3 ?
+ (i<4*fe.dofs_per_vertex ?
+ i
:
- (dim == 3 ?
- (i<4*fe.dofs_per_vertex ?
- i
- :
- (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
- i+4*fe.dofs_per_vertex
- :
- i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+ (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
+ i+4*fe.dofs_per_vertex
:
- numbers::invalid_unsigned_int)));
+ i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+ :
+ numbers::invalid_unsigned_int)));
if (fe.is_primitive (cell_index))
- {
+ {
if (component_mask[fe.face_system_to_component_index(i).first]
== true)
selected_dofs.add_index (face_dof_indices[i]);
continue;
for (unsigned int dof=0; dof!=fe.dofs_per_line; ++dof)
dofs_on_children.push_back(
- this_face->child(child)->dof_index(dof, fe_index));
+ this_face->child(child)->dof_index(dof, fe_index));
}
// note: can get fewer DoFs when we have artificial cells
Assert(dofs_on_children.size() <= n_dofs_on_children, ExcInternalError());
if (cell->neighbor_child_on_subface (face, child)->is_artificial())
continue;
for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof)
- dofs_on_children.push_back(
- this_face->child(child)->dof_index(dof, fe_index));
+ dofs_on_children.push_back(
+ this_face->child(child)->dof_index(dof, fe_index));
}
// note: can get fewer DoFs when we have artificial cells:
// there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs.
// FESystem(FE_Q(1),FE_DGQ(1)). neither of them
// dominates the other.
- //
- // a final possibility is that we have something like
- // FESystem(FE_Q(1),FE_Q(1)) vs
- // FESystem(FE_Q(1),FE_Nothing()), see
- // hp/fe_nothing_18/19.
- //
- // in any case, the point is that it doesn't
+ //
+ // a final possibility is that we have something like
+ // FESystem(FE_Q(1),FE_Q(1)) vs
+ // FESystem(FE_Q(1),FE_Nothing()), see
+ // hp/fe_nothing_18/19.
+ //
+ // in any case, the point is that it doesn't
// matter. there is nothing to do here.
break;
}
face_1->get_dof_indices(dofs_1, face_1_index);
face_2->get_dof_indices(dofs_2, face_2_index);
- for (unsigned int i=0; i < dofs_per_face; i++) {
- if (dofs_1[i] == numbers::invalid_dof_index ||
- dofs_2[i] == numbers::invalid_dof_index) {
- /* If either of these faces have no indices, stop. This is so
- * that there is no attempt to match artificial cells of
- * parallel distributed triangulations.
- *
- * While it seems like we ought to be able to avoid even calling
- * set_periodicity_constraints for artificial faces, this
- * situation can arise when a face that is being made periodic
- * is only partially touched by the local subdomain.
- * make_periodicity_constraints will be called recursively even
- * for the section of the face that is not touched by the local
- * subdomain.
- *
- * Until there is a better way to determine if the cells that
- * neighbor a face are artificial, we simply test to see if the
- * face does not have a valid dof initialization.
- */
- return;
+ for (unsigned int i=0; i < dofs_per_face; i++)
+ {
+ if (dofs_1[i] == numbers::invalid_dof_index ||
+ dofs_2[i] == numbers::invalid_dof_index)
+ {
+ /* If either of these faces have no indices, stop. This is so
+ * that there is no attempt to match artificial cells of
+ * parallel distributed triangulations.
+ *
+ * While it seems like we ought to be able to avoid even calling
+ * set_periodicity_constraints for artificial faces, this
+ * situation can arise when a face that is being made periodic
+ * is only partially touched by the local subdomain.
+ * make_periodicity_constraints will be called recursively even
+ * for the section of the face that is not touched by the local
+ * subdomain.
+ *
+ * Until there is a better way to determine if the cells that
+ * neighbor a face are artificial, we simply test to see if the
+ * face does not have a valid dof initialization.
+ */
+ return;
+ }
}
- }
// Well, this is a hack:
//
// orientation:
const unsigned int j =
cell_to_rotated_face_index[fe.face_to_cell_index(identity_constraint_target,
- 0, /* It doesn't really matter, just assume
+ 0, /* It doesn't really matter, just assume
* we're on the first face...
*/
face_orientation, face_flip, face_rotation)];
"different to denote different boundaries."));
std::vector<GridTools::PeriodicFacePair
- <typename DH::cell_iterator> > matched_faces;
+ <typename DH::cell_iterator> > matched_faces;
// Collect matching periodic cells on the coarsest level:
- GridTools::collect_periodic_faces(dof_handler, b_id1, b_id2, direction,
+ GridTools::collect_periodic_faces(dof_handler, b_id1, b_id2, direction,
matched_faces, offset);
make_periodicity_constraints<DH>
- (matched_faces, constraint_matrix, component_mask);
+ (matched_faces, constraint_matrix, component_mask);
}
ExcNotImplemented());
std::vector<GridTools::PeriodicFacePair
- <typename DH::cell_iterator> > matched_faces;
+ <typename DH::cell_iterator> > matched_faces;
// Collect matching periodic cells on the coarsest level:
GridTools::collect_periodic_faces(dof_handler, b_id, direction,
matched_faces, offset);
make_periodicity_constraints<DH>
- (matched_faces, constraint_matrix, component_mask);
+ (matched_faces, constraint_matrix, component_mask);
}
void
make_periodicity_constraints
(const std::vector<GridTools::PeriodicFacePair<typename DH::cell_iterator> >
- &periodic_faces,
+ &periodic_faces,
dealii::ConstraintMatrix &constraint_matrix,
const ComponentMask &component_mask)
{
// And apply the low level make_periodicity_constraints function to
// every matching pair:
- for(; it!=end_periodic; ++it)
- {
- typedef typename DH::face_iterator FaceIterator;
- const FaceIterator face_1 = it->cell[0]->face(it->face_idx[0]);
- const FaceIterator face_2 = it->cell[1]->face(it->face_idx[1]);
+ for (; it!=end_periodic; ++it)
+ {
+ typedef typename DH::face_iterator FaceIterator;
+ const FaceIterator face_1 = it->cell[0]->face(it->face_idx[0]);
+ const FaceIterator face_2 = it->cell[1]->face(it->face_idx[1]);
- Assert(face_1->at_boundary() && face_2->at_boundary(),
- ExcInternalError());
+ Assert(face_1->at_boundary() && face_2->at_boundary(),
+ ExcInternalError());
- Assert (face_1 != face_2,
- ExcInternalError());
+ Assert (face_1 != face_2,
+ ExcInternalError());
- make_periodicity_constraints(face_1,
- face_2,
- constraint_matrix,
- component_mask,
- it->orientation[0],
- it->orientation[1],
- it->orientation[2]);
- }
+ make_periodicity_constraints(face_1,
+ face_2,
+ constraint_matrix,
+ component_mask,
+ it->orientation[0],
+ it->orientation[1],
+ it->orientation[2]);
+ }
}
*/
template <int dim, int spacedim>
void compute_intergrid_weights_3 (
- const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &cell,
- const Assembler::Scratch &,
- Assembler::CopyData<dim,spacedim> ©_data,
- const unsigned int coarse_component,
- const FiniteElement<dim,spacedim> &coarse_fe,
- const InterGridMap<dealii::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
- const std::vector<dealii::Vector<double> > ¶meter_dofs,
- const std::vector<types::global_dof_index> &weight_mapping)
+ const typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator &cell,
+ const Assembler::Scratch &,
+ Assembler::CopyData<dim,spacedim> ©_data,
+ const unsigned int coarse_component,
+ const FiniteElement<dim,spacedim> &coarse_fe,
+ const InterGridMap<dealii::DoFHandler<dim,spacedim> > &coarse_to_fine_grid_map,
+ const std::vector<dealii::Vector<double> > ¶meter_dofs,
+ const std::vector<types::global_dof_index> &weight_mapping)
{
// for each cell on the parameter grid: find out which degrees of
// freedom on the fine grid correspond in which way to the degrees
// parameter grid cell
cell->get_dof_indices (copy_data.parameter_dof_indices);
- // reset the output array to a pristine state
- copy_data.global_parameter_representation.clear ();
+ // reset the output array to a pristine state
+ copy_data.global_parameter_representation.clear ();
// loop over all dofs on this cell and check whether they are
// interesting for us
= coarse_fe.system_to_component_index(local_dof).second;
copy_data.global_parameter_representation.push_back(
- dealii::Vector<double> (n_fine_dofs));
+ dealii::Vector<double> (n_fine_dofs));
// distribute the representation of
// @p{local_parameter_dof} on the parameter grid cell
*/
template <int dim,int spacedim>
void copy_intergrid_weights_3(const Assembler::CopyData<dim,spacedim> ©_data,
- const unsigned int coarse_component,
+ const unsigned int coarse_component,
const FiniteElement<dim,spacedim> &coarse_fe,
- const std::vector<types::global_dof_index> &weight_mapping,
+ const std::vector<types::global_dof_index> &weight_mapping,
std::vector<std::map<types::global_dof_index, float> > &weights)
{
unsigned int pos = 0;
if (coarse_fe.system_to_component_index(local_dof).first
==
coarse_component)
- {
+ {
// now that we've got the global representation of each
// parameter dof, we've only got to clobber the non-zero
// entries in that vector and store the result
// several threads operating on different intergrid
// weights, have only one mutex for all of them
for (types::global_dof_index i=0; i<copy_data.global_parameter_representation[pos].size();
- ++i)
+ ++i)
// set this weight if it belongs to a parameter dof.
if (weight_mapping[i] != numbers::invalid_dof_index)
{
Assembler::CopyData<dim,spacedim> copy_data;
WorkStream::run(coarse_grid.begin_active(),
- coarse_grid.end(),
- std_cxx1x::bind(&compute_intergrid_weights_3<dim,spacedim>,
- std_cxx1x::_1,
- std_cxx1x::_2,
- std_cxx1x::_3,
- coarse_component,
- std_cxx1x::cref(coarse_grid.get_fe()),
- std_cxx1x::cref(coarse_to_fine_grid_map),
- std_cxx1x::cref(parameter_dofs),
- std_cxx1x::cref(weight_mapping)),
- std_cxx1x::bind(©_intergrid_weights_3<dim,spacedim>,
- std_cxx1x::_1,
- coarse_component,
- std_cxx1x::cref(coarse_grid.get_fe()),
- std_cxx1x::cref(weight_mapping),
- std_cxx1x::ref(weights)),
- scratch,
- copy_data);
+ coarse_grid.end(),
+ std_cxx1x::bind(&compute_intergrid_weights_3<dim,spacedim>,
+ std_cxx1x::_1,
+ std_cxx1x::_2,
+ std_cxx1x::_3,
+ coarse_component,
+ std_cxx1x::cref(coarse_grid.get_fe()),
+ std_cxx1x::cref(coarse_to_fine_grid_map),
+ std_cxx1x::cref(parameter_dofs),
+ std_cxx1x::cref(weight_mapping)),
+ std_cxx1x::bind(©_intergrid_weights_3<dim,spacedim>,
+ std_cxx1x::_1,
+ coarse_component,
+ std_cxx1x::cref(coarse_grid.get_fe()),
+ std_cxx1x::cref(weight_mapping),
+ std_cxx1x::ref(weights)),
+ scratch,
+ copy_data);
}
if (!cell->is_artificial()
&&
cell->at_boundary ())
- {
- const FiniteElement<dim,spacedim> &fe = cell->get_fe();
+ {
+ const FiniteElement<dim,spacedim> &fe = cell->get_fe();
- // get global indices of dofs on the cell
- cell_dofs.resize (fe.dofs_per_cell);
- cell->get_dof_indices (cell_dofs);
+ // get global indices of dofs on the cell
+ cell_dofs.resize (fe.dofs_per_cell);
+ cell->get_dof_indices (cell_dofs);
- for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
- ++face_no)
- {
- const typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
-
- // if face is on the boundary and satisfies the correct
- // boundary id property
- if (face->at_boundary ()
- &&
- ((boundary_indicator == numbers::invalid_boundary_id)
- ||
- (face->boundary_indicator() == boundary_indicator)))
- {
- // get indices and physical location on this face
- face_dofs.resize (fe.dofs_per_face);
- face->get_dof_indices (face_dofs, cell->active_fe_index());
+ for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
+ ++face_no)
+ {
+ const typename DH<dim,spacedim>::face_iterator face = cell->face(face_no);
+
+ // if face is on the boundary and satisfies the correct
+ // boundary id property
+ if (face->at_boundary ()
+ &&
+ ((boundary_indicator == numbers::invalid_boundary_id)
+ ||
+ (face->boundary_indicator() == boundary_indicator)))
+ {
+ // get indices and physical location on this face
+ face_dofs.resize (fe.dofs_per_face);
+ face->get_dof_indices (face_dofs, cell->active_fe_index());
- // enter those dofs into the list that match the component
- // signature.
- for (unsigned int i=0; i<face_dofs.size(); ++i)
- {
- // Find out if a dof has a contribution in this
- // component, and if so, add it to the list
- const std::vector<types::global_dof_index>::iterator it_index_on_cell
- = std::find (cell_dofs.begin(), cell_dofs.end(), face_dofs[i]);
- Assert (it_index_on_cell != cell_dofs.end(), ExcInvalidIterator());
- const unsigned int index_on_cell = std::distance(cell_dofs.begin(),
- it_index_on_cell);
- const ComponentMask &nonzero_component_array
- = cell->get_fe().get_nonzero_components (index_on_cell);
- bool nonzero = false;
- for (unsigned int c=0; c<n_components; ++c)
- if (nonzero_component_array[c] && component_mask[c])
- {
- nonzero = true;
- break;
- }
+ // enter those dofs into the list that match the component
+ // signature.
+ for (unsigned int i=0; i<face_dofs.size(); ++i)
+ {
+ // Find out if a dof has a contribution in this
+ // component, and if so, add it to the list
+ const std::vector<types::global_dof_index>::iterator it_index_on_cell
+ = std::find (cell_dofs.begin(), cell_dofs.end(), face_dofs[i]);
+ Assert (it_index_on_cell != cell_dofs.end(), ExcInvalidIterator());
+ const unsigned int index_on_cell = std::distance(cell_dofs.begin(),
+ it_index_on_cell);
+ const ComponentMask &nonzero_component_array
+ = cell->get_fe().get_nonzero_components (index_on_cell);
+ bool nonzero = false;
+ for (unsigned int c=0; c<n_components; ++c)
+ if (nonzero_component_array[c] && component_mask[c])
+ {
+ nonzero = true;
+ break;
+ }
- if (nonzero)
- zero_boundary_constraints.add_line (face_dofs[i]);
- }
- }
- }
- }
+ if (nonzero)
+ zero_boundary_constraints.add_line (face_dofs[i]);
+ }
+ }
+ }
+ }
}
system_to_component_table.resize(this->dofs_per_cell);
face_system_to_component_table.resize(this->dofs_per_face);
for (unsigned int j=0 ; j<this->dofs_per_cell ; ++j)
- system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
+ system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
for (unsigned int j=0 ; j<this->dofs_per_face ; ++j)
- face_system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
+ face_system_to_component_table[j] = std::pair<unsigned,unsigned>(0,j);
}
-
+
for (unsigned int j=0 ; j<this->dofs_per_cell ; ++j)
system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);
for (unsigned int j=0 ; j<this->dofs_per_face ; ++j)
face_system_to_base_table[j] = std::make_pair(std::make_pair(0U,0U),j);
-
+
// Fill with default value; may be changed by constructor of derived class.
base_to_block_indices.reinit(1,1);
{
switch (direction)
{
- // Rotate xy-plane
- // counter-clockwise
+ // Rotate xy-plane
+ // counter-clockwise
case 'z':
for (unsigned int iz=0; iz<((dim>2) ? n:1); ++iz)
for (unsigned int j=0; j<n; ++j)
numbers[l++] = k;
}
break;
- // Rotate xy-plane
- // clockwise
+ // Rotate xy-plane
+ // clockwise
case 'Z':
for (unsigned int iz=0; iz<((dim>2) ? n:1); ++iz)
for (unsigned int iy=0; iy<n; ++iy)
numbers[k] = l++;
}
break;
- // Rotate yz-plane
- // counter-clockwise
+ // Rotate yz-plane
+ // counter-clockwise
case 'x':
Assert (dim>2, ExcDimensionMismatch (dim,3));
for (unsigned int iz=0; iz<n; ++iz)
numbers[l++] = k;
}
break;
- // Rotate yz-plane
- // clockwise
+ // Rotate yz-plane
+ // clockwise
case 'X':
Assert (dim>2, ExcDimensionMismatch (dim,3));
for (unsigned int iz=0; iz<n; ++iz)
for (unsigned int iy=0; iy <= ((codim>1) ? this->degree : 0) ; ++iy)
for (unsigned int ix=0; ix<=this->degree; ++ix)
{
- Point<codim> p;
-
+ Point<codim> p;
+
p(0) = points[ix][0];
if (codim>1)
p(1) = points[iy][0];
(typename FiniteElement<dim,spacedim>::
ExcInterpolationNotImplemented ()));
- // generate a quadrature with the unit face support points.
+ // generate a quadrature with the unit face support points.
const Quadrature<dim-1> face_quadrature (source_fe->get_unit_face_support_points ());
// Rule of thumb for FP accuracy, that can be expected for a given
:
FE_PolyFace<PolynomialSpace<dim-1>, dim, spacedim>
(PolynomialSpace<dim-1>(Polynomials::Legendre::generate_complete_basis(degree)),
- FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
- std::vector<bool>(1,true))
+ FiniteElementData<dim>(get_dpo_vector(degree), 1, degree, FiniteElementData<dim>::L2),
+ std::vector<bool>(1,true))
{}
const double eps = 2e-13*(this->degree+1)*(dim-1);
FullMatrix<double> mass (face_quadrature.size(), source_fe->dofs_per_face);
-
+
for (unsigned int k = 0; k < face_quadrature.size(); ++k)
{
const Point<dim-1> p =
for (unsigned int k = 0; k < face_quadrature.size(); ++k)
{
const Point<dim-1> p = numbers::invalid_unsigned_int ?
- face_quadrature.point(k) :
- GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(k),
- subface);
+ face_quadrature.point(k) :
+ GeometryInfo<dim-1>::child_to_cell_coordinates (face_quadrature.point(k),
+ subface);
v_in(k) = this->poly_space.compute_value(i, p);
}
const double result = H.least_squares(v_out, v_in);
// be done for the face orientation flag in 3D.
const Quadrature<dim> subface_quadrature
= subface == numbers::invalid_unsigned_int
- ?
- QProjector<dim>::project_to_face (quad_face_support, 0)
- :
- QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
+ ?
+ QProjector<dim>::project_to_face (quad_face_support, 0)
+ :
+ QProjector<dim>::project_to_subface (quad_face_support, 0, subface);
for (unsigned int i=0; i<source_fe->dofs_per_face; ++i)
{
const Point<dim> &p = subface_quadrature.point (i);
// face orientation, face flips, etc
unsigned int adjusted_dof_index_on_line;
switch (dim)
- {
- case 1:
- Assert (false, ExcInternalError());
-
- case 2:
- // in 2d, only face_flip has a meaning. if it is set, consider
- // dofs in reverse order
- if (face_flip == false)
+ {
+ case 1:
+ Assert (false, ExcInternalError());
+
+ case 2:
+ // in 2d, only face_flip has a meaning. if it is set, consider
+ // dofs in reverse order
+ if (face_flip == false)
+ adjusted_dof_index_on_line = dof_index_on_line;
+ else
+ adjusted_dof_index_on_line = this->dofs_per_line - dof_index_on_line - 1;
+ break;
+
+ case 3:
+ // in 3d, things are difficult. someone will have to think
+ // about how this code here should look like, by drawing a bunch
+ // of pictures of how all the faces can look like with the various
+ // flips and rotations.
+ //
+ // that said, the Q2 case is easy enough to implement, as is the case
+ // where everything is in standard orientation
+ Assert ((this->dofs_per_line <= 1) ||
+ ((face_orientation == true) &&
+ (face_flip == false) &&
+ (face_rotation == false)),
+ ExcNotImplemented());
adjusted_dof_index_on_line = dof_index_on_line;
- else
- adjusted_dof_index_on_line = this->dofs_per_line - dof_index_on_line - 1;
- break;
-
- case 3:
- // in 3d, things are difficult. someone will have to think
- // about how this code here should look like, by drawing a bunch
- // of pictures of how all the faces can look like with the various
- // flips and rotations.
- //
- // that said, the Q2 case is easy enough to implement, as is the case
- // where everything is in standard orientation
- Assert ((this->dofs_per_line <= 1) ||
- ((face_orientation == true) &&
- (face_flip == false) &&
- (face_rotation == false)),
- ExcNotImplemented());
- adjusted_dof_index_on_line = dof_index_on_line;
- break;
- }
+ break;
+ }
return (this->first_line_index
+ GeometryInfo<dim>::face_to_cell_lines(face, face_line,
template <int dim, int spacedim>
FESystem<dim,spacedim>::InternalData::InternalData(const unsigned int n_base_elements,
const bool compute_hessians)
-:
+ :
compute_hessians (compute_hessians),
base_fe_datas(n_base_elements),
base_fe_values_datas(n_base_elements)
const unsigned int
base_face_to_cell_index
- = this->base_element(face_base_index.first.first).face_to_cell_index (face_base_index.second,
- face,
- face_orientation,
- face_flip,
- face_rotation);
+ = this->base_element(face_base_index.first.first).face_to_cell_index (face_base_index.second,
+ face,
+ face_orientation,
+ face_flip,
+ face_rotation);
// it would be nice if we had a base_to_system_index function, but
// all that exists is a component_to_system_index function. we can't do
start_indices[i] = k;
k += block_data[i];
}
-
+
for (unsigned int i=0; i<element.dofs_per_cell; ++i)
{
std::pair<unsigned int, types::global_dof_index>
namespace FETools
{
template <int dim, int spacedim,
- template <int, int> class DH1,
- template <int, int> class DH2,
- class InVector, class OutVector>
+ template <int, int> class DH1,
+ template <int, int> class DH2,
+ class InVector, class OutVector>
void
interpolate(const DH1<dim, spacedim> &dof1,
const InVector &u1,
template <int dim, int spacedim,
- template <int, int> class DH1,
- template <int, int> class DH2,
- class InVector, class OutVector>
+ template <int, int> class DH1,
+ template <int, int> class DH2,
+ class InVector, class OutVector>
void
interpolate (const DH1<dim, spacedim> &dof1,
const InVector &u1,
template <int dim,
- template <int> class DH,
- class InVector, class OutVector, int spacedim>
+ template <int> class DH,
+ class InVector, class OutVector, int spacedim>
void
back_interpolate(const DH<dim> &dof1,
const InVector &u1,
(typeid(*present_cell.get()) == typeid(Type)))
{
// call destructor of the old object
- static_cast<const Type*>(present_cell.get())->~Type();
+ static_cast<const Type *>(present_cell.get())->~Type();
// then construct a new object in-place
new(const_cast<void *>(static_cast<const void *>(present_cell.get()))) Type(new_cell);
// Qp mapping of cells at the boundary.
const QGaussLobatto<1> points(degree+1);
tensor_pols = new TensorProductPolynomials<dim>
- (Polynomials::generate_complete_Lagrange_basis(points.get_points()));
+ (Polynomials::generate_complete_Lagrange_basis(points.get_points()));
Assert (n_shape_functions==tensor_pols->n(),
ExcInternalError());
Assert(n_inner+n_outer==n_shape_functions, ExcInternalError());
double const *loqv_ptr=0;
switch (degree)
{
- // for degree==1, we shouldn't have to compute any support points, since
- // all of them are on the vertices
+ // for degree==1, we shouldn't have to compute any support points, since
+ // all of them are on the vertices
case 2:
{
for (unsigned int point=0; point<n_q_points; ++point)
for (unsigned int i=0; i<n_inner; ++i)
for (unsigned int j=0; j<n_inner; ++j)
- {
- long double res = 0.;
- for (unsigned int l=0; l<dim; ++l)
- res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
- (long double)quadrature_data.derivative(point, n_outer+j)[l];
+ {
+ long double res = 0.;
+ for (unsigned int l=0; l<dim; ++l)
+ res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
+ (long double)quadrature_data.derivative(point, n_outer+j)[l];
- S(i,j) += res * (long double)quadrature.weight(point);
- }
+ S(i,j) += res * (long double)quadrature.weight(point);
+ }
// Compute the components of T to be the product of gradients of inner and
// outer shape functions.
for (unsigned int point=0; point<n_q_points; ++point)
for (unsigned int i=0; i<n_inner; ++i)
for (unsigned int k=0; k<n_outer; ++k)
- {
- long double res = 0.;
- for (unsigned int l=0; l<dim; ++l)
- res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
- (long double)quadrature_data.derivative(point, k)[l];
-
- T(i,k) += res *(long double)quadrature.weight(point);
- }
+ {
+ long double res = 0.;
+ for (unsigned int l=0; l<dim; ++l)
+ res += (long double)quadrature_data.derivative(point, n_outer+i)[l] *
+ (long double)quadrature_data.derivative(point, k)[l];
+
+ T(i,k) += res *(long double)quadrature.weight(point);
+ }
FullMatrix<long double> S_1(n_inner);
S_1.invert(S);
// early out if we already have our point
if (f.square() < 1e-24 * cell->diameter() * cell->diameter())
return p_unit;
-
+
// we need to compare the position of the computed p(x) against the given
// point 'p'. We will terminate the iteration and return 'x' if they are
// less than eps apart. The question is how to choose eps -- or, put maybe
#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
std::cout << "Newton iteration " << newton_iteration << std::endl;
#endif
-
+
// f'(x)
Tensor<2,spacedim> df;
for (unsigned int k=0; k<mdata.n_shape_functions; ++k)
#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
std::cout << " delta=" << delta << std::endl;
#endif
-
+
// do a line search
double step_length = 1;
do
const Point<spacedim> f_trial = p_real_trial-p;
#ifdef DEBUG_TRANSFORM_REAL_TO_UNIT_CELL
- std::cout << " step_length=" << step_length << std::endl
- << " ||f || =" << f.norm() << std::endl
- << " ||f*|| =" << f_trial.norm() << std::endl
- << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl;
+ std::cout << " step_length=" << step_length << std::endl
+ << " ||f || =" << f.norm() << std::endl
+ << " ||f*|| =" << f_trial.norm() << std::endl
+ << " ||f*||_A =" << (df_inverse * f_trial).norm() << std::endl;
#endif
-
+
// see if we are making progress with the current step length
// and if not, reduce it by a factor of two and try again
- //
- // strictly speaking, we should probably use the same norm as we use
- // for the outer algorithm. in practice, line search is just a
- // crutch to find a "reasonable" step length, and so using the l2
- // norm is probably just fine
+ //
+ // strictly speaking, we should probably use the same norm as we use
+ // for the outer algorithm. in practice, line search is just a
+ // crutch to find a "reasonable" step length, and so using the l2
+ // norm is probably just fine
if (f_trial.norm() < f.norm())
{
p_real = p_real_trial;
else if (step_length > 0.05)
step_length /= 2;
else
- AssertThrow (false,
- (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
+ AssertThrow (false,
+ (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
}
while (true);
++newton_iteration;
if (newton_iteration > newton_iteration_limit)
- AssertThrow (false,
- (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
+ AssertThrow (false,
+ (typename Mapping<dim,spacedim>::ExcTransformationFailed()));
last_f_weighted_norm = (df_inverse * f).norm();
}
while (last_f_weighted_norm > eps);
-
+
return p_unit;
}
switch (dim)
{
- // A line (1d parallelepiped)
+ // A line (1d parallelepiped)
case 1:
vertices[1] = corners[0];
break;
- // A parallelogram (2d parallelepiped)
+ // A parallelogram (2d parallelepiped)
case 2:
// assign corners to vertices:
vertices[1] = corners[0];
vertices[3] = vertices[1] + vertices[2];
break;
- // A parallelepiped (3d parallelepiped)
+ // A parallelepiped (3d parallelepiped)
case 3:
// assign corners to vertices:
vertices[1] = corners[0];
// Determine number of cells and vertices
const unsigned int
n_cells = static_cast<unsigned int>(std::ceil (half_length /
- std::max (radius_0,
- radius_1)));
+ std::max (radius_0,
+ radius_1)));
const unsigned int n_vertices = 4 * (n_cells + 1);
std::vector<Point<3> > vertices_tmp(n_vertices);
text[2] = "ASCII";
text[3] = "DATASET UNSTRUCTURED_GRID";
- for(unsigned int i = 0;i < 4;i++)//Checking for the match between initial strings/text in the file.
+ for (unsigned int i = 0; i < 4; i++) //Checking for the match between initial strings/text in the file.
{
getline(in,line);
in>>no_vertices;// taking the no. of vertices
in.ignore(256, '\n');//ignoring the number beside the total no. of points.
- for(unsigned int count = 0; count < no_vertices; count++)//loop to read three values till the no . vertices is satisfied
+ for (unsigned int count = 0; count < no_vertices; count++) //loop to read three values till the no . vertices is satisfied
{
- // VTK format always specifies vertex coordinates with 3 components
- Point<3> x;
+ // VTK format always specifies vertex coordinates with 3 components
+ Point<3> x;
in >> x(0) >> x(1) >> x(2);
vertices.push_back(Point<spacedim>());
- for (unsigned int d=0; d<spacedim; ++d)
- vertices.back()(d) = x(d);
+ for (unsigned int d=0; d<spacedim; ++d)
+ vertices.back()(d) = x(d);
vertex_indices[count] = count;
}
if (dim == 3)
{
- for(unsigned int count = 0; count < total_cells; count++)
+ for (unsigned int count = 0; count < total_cells; count++)
{
in>>type;
- if(type == 8)
+ if (type == 8)
{
cells.push_back(CellData<dim>());
- for(j = 0; j < type; j++)//loop to feed data
+ for (j = 0; j < type; j++) //loop to feed data
in >> cells.back().vertices[j];
cells.back().material_id = 0;
- for(j = 0; j < type; j++)//loop to feed the data of the vertices to the cell
+ for (j = 0; j < type; j++) //loop to feed the data of the vertices to the cell
{
cells.back().vertices[j] = vertex_indices[cells.back().vertices[j]];
}
subcelldata.boundary_quads.push_back(CellData<2>());
int no_quads = 0;
- for(j = 0; j < type; j++)//loop to feed the data to the boundary
+ for (j = 0; j < type; j++) //loop to feed the data to the boundary
{
in >> subcelldata.boundary_quads.back().vertices[j];
}
subcelldata.boundary_quads.back().material_id = 0;
- for(unsigned int j = 0; j < type; j++)
+ for (unsigned int j = 0; j < type; j++)
{
subcelldata.boundary_quads.back().vertices[j] = vertex_indices[subcelldata.boundary_quads.back().vertices[j]];
}
else if (dim == 2)
{
- for(unsigned int count = 0; count < total_cells; count++)
+ for (unsigned int count = 0; count < total_cells; count++)
{
in>>type;
- if(type == 4)
+ if (type == 4)
{
cells.push_back(CellData<dim>());
- for(j = 0; j < type; j++)//loop to feed data
+ for (j = 0; j < type; j++) //loop to feed data
in >> cells.back().vertices[j];
cells.back().material_id = 0;
- for(j = 0; j < type; j++)//loop to feed the data of the vertices to the cell
+ for (j = 0; j < type; j++) //loop to feed the data of the vertices to the cell
{
cells.back().vertices[j] = vertex_indices[cells.back().vertices[j]];
}
no_cells++;
}
- else if(type == 2)
+ else if (type == 2)
{
//If this is encountered, the pointer comes out of the loop
//and starts processing boundaries.
subcelldata.boundary_lines.push_back(CellData<1>());
int no_lines = 0;
- for(j = 0; j < type; j++)//loop to feed the data to the boundary
+ for (j = 0; j < type; j++) //loop to feed the data to the boundary
{
in >> subcelldata.boundary_lines.back().vertices[j];
}
subcelldata.boundary_lines.back().material_id = 0;
- for(unsigned int j = 0; j < type; j++)
+ for (unsigned int j = 0; j < type; j++)
{
subcelldata.boundary_lines.back().vertices[j] = vertex_indices[subcelldata.boundary_lines.back().vertices[j]];
}
{
in.ignore(256, '\n');
- while(!in.eof())
+ while (!in.eof())
{
in>>keyword_3;
- if(std::strcmp(keyword_3,"12") && std::strcmp(keyword_3,"9"))
+ if (std::strcmp(keyword_3,"12") && std::strcmp(keyword_3,"9"))
{
break;
}
in.ignore(256, '\n');
- for(unsigned int i = 0;i < 2;i++)
+ for (unsigned int i = 0; i < 2; i++)
{
getline(in,linenew);
- if(linenew.compare(textnew[i]) == 0)
+ if (linenew.compare(textnew[i]) == 0)
{}
else
AssertThrow (false,
textnew[i] + "> section"));
}
- for(unsigned int i = 0; i < no_cells; i++)//assigning IDs to cells.
+ for (unsigned int i = 0; i < no_cells; i++) //assigning IDs to cells.
{
int id;
in>>id;
if (dim == 3)
{
- for(unsigned int i = 0; i < no_quads; i++)//assigning IDs to bounds.
+ for (unsigned int i = 0; i < no_quads; i++) //assigning IDs to bounds.
{
int id;
in>>id;
}
else if (dim == 2)
{
- for(unsigned int i = 0; i < no_lines; i++)//assigning IDs to bounds.
+ for (unsigned int i = 0; i < no_lines; i++) //assigning IDs to bounds.
{
int id;
in>>id;
cells,
subcelldata);
- if(dim == spacedim)
+ if (dim == spacedim)
GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(vertices,
- cells);
+ cells);
GridReordering<dim, spacedim>::reorder_cells(cells);
tria->create_triangulation_compatibility(vertices,
return;
case vtk:
- read_vtk (in);
- return;
+ read_vtk (in);
+ return;
case unv:
read_unv (in);
for (unsigned int i=0; i<vertices.size(); ++i)
if (vertex_used[i])
- out << '\t' << vertices[i] << '\n';
+ out << '\t' << vertices[i] << '\n';
// write cells or faces
const bool write_cells = dx_flags.write_cells;
}
AssertThrow (best_cell.first.state() == IteratorState::valid,
- ExcPointNotFound<spacedim>(p));
+ ExcPointNotFound<spacedim>(p));
return best_cell;
}
}
AssertThrow (best_cell.first.state() == IteratorState::valid,
- ExcPointNotFound<spacedim>(p));
+ ExcPointNotFound<spacedim>(p));
return best_cell;
}
// pointer/index clashes when saving/restoring them. The following approach
// works, but this map can get quite big. Not sure about more efficient solutions.
std::map< std::pair<unsigned int,unsigned int>, unsigned int >
- indexmap;
+ indexmap;
unsigned int index = 0;
for (typename Triangulation<dim,spacedim>::active_cell_iterator
cell = triangulation.begin_active();
(cell->neighbor(f)->has_children() == false))
{
unsigned int other_index = indexmap.find(
- std::pair<unsigned int,unsigned int>(cell->neighbor(f)->level(),cell->neighbor(f)->index()))->second;
+ std::pair<unsigned int,unsigned int>(cell->neighbor(f)->level(),cell->neighbor(f)->index()))->second;
cell_connectivity.add (index, other_index);
cell_connectivity.add (other_index, index);
}
*/
template<typename CellIterator>
void
- match_periodic_face_pairs
- (std::set<std::pair<CellIterator, unsigned int> > &pairs1,
- std::set<std::pair<typename identity<CellIterator>::type, unsigned int> > &pairs2,
- const int direction,
- std::vector<PeriodicFacePair<CellIterator> > & matched_pairs,
- const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset)
+ match_periodic_face_pairs
+ (std::set<std::pair<CellIterator, unsigned int> > &pairs1,
+ std::set<std::pair<typename identity<CellIterator>::type, unsigned int> > &pairs2,
+ const int direction,
+ std::vector<PeriodicFacePair<CellIterator> > &matched_pairs,
+ const dealii::Tensor<1,CellIterator::AccessorType::space_dimension> &offset)
{
static const int space_dim = CellIterator::AccessorType::space_dimension;
Assert (0<=direction && direction<space_dim,
// Match with a complexity of O(n^2). This could be improved...
std::bitset<3> orientation;
typedef typename std::set
- <std::pair<CellIterator, unsigned int> >::const_iterator PairIterator;
+ <std::pair<CellIterator, unsigned int> >::const_iterator PairIterator;
for (PairIterator it1 = pairs1.begin(); it1 != pairs1.end(); ++it1)
{
for (PairIterator it2 = pairs2.begin(); it2 != pairs2.end(); ++it2)
// remove the matched cell in pairs2 to speed up the
// matching:
const PeriodicFacePair<CellIterator> matched_face
- = {{cell1, cell2},{face_idx1, face_idx2}, orientation};
+ = {{cell1, cell2},{face_idx1, face_idx2}, orientation};
matched_pairs.push_back(matched_face);
pairs2.erase(it2);
++n_matches;
template<typename CONTAINER>
void
collect_periodic_faces
- (const CONTAINER &container,
- const types::boundary_id b_id1,
- const types::boundary_id b_id2,
- const int direction,
- std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > & matched_pairs,
- const dealii::Tensor<1,CONTAINER::space_dimension> &offset)
+ (const CONTAINER &container,
+ const types::boundary_id b_id1,
+ const types::boundary_id b_id2,
+ const int direction,
+ std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
+ const dealii::Tensor<1,CONTAINER::space_dimension> &offset)
{
static const int dim = CONTAINER::dimension;
static const int space_dim = CONTAINER::space_dimension;
{
const typename CONTAINER::face_iterator face = cell->face(i);
if (face->at_boundary() && face->boundary_indicator() == b_id1)
- {
- const std::pair<typename CONTAINER::cell_iterator, unsigned int> pair1
- = std::make_pair(cell, i);
- pairs1.insert(pair1);
- }
+ {
+ const std::pair<typename CONTAINER::cell_iterator, unsigned int> pair1
+ = std::make_pair(cell, i);
+ pairs1.insert(pair1);
+ }
if (face->at_boundary() && face->boundary_indicator() == b_id2)
- {
- const std::pair<typename CONTAINER::cell_iterator, unsigned int> pair2
- = std::make_pair(cell, i);
- pairs2.insert(pair2);
- }
+ {
+ const std::pair<typename CONTAINER::cell_iterator, unsigned int> pair2
+ = std::make_pair(cell, i);
+ pairs2.insert(pair2);
+ }
}
}
template<typename CONTAINER>
void
collect_periodic_faces
- (const CONTAINER &container,
- const types::boundary_id b_id,
- const int direction,
- std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > & matched_pairs,
- const dealii::Tensor<1,CONTAINER::space_dimension> &offset)
+ (const CONTAINER &container,
+ const types::boundary_id b_id,
+ const int direction,
+ std::vector<PeriodicFacePair<typename CONTAINER::cell_iterator> > &matched_pairs,
+ const dealii::Tensor<1,CONTAINER::space_dimension> &offset)
{
static const int dim = CONTAINER::dimension;
static const int space_dim = CONTAINER::space_dimension;
for (typename CONTAINER::cell_iterator cell = container.begin(0);
cell != container.end(0); ++cell)
- {
- const typename CONTAINER::face_iterator face_1 = cell->face(2*direction);
- const typename CONTAINER::face_iterator face_2 = cell->face(2*direction+1);
-
- if (face_1->at_boundary() && face_1->boundary_indicator() == b_id)
- {
- const std::pair<typename CONTAINER::cell_iterator, unsigned int> pair1
- = std::make_pair(cell, 2*direction);
- pairs1.insert(pair1);
- }
-
- if (face_2->at_boundary() && face_2->boundary_indicator() == b_id)
- {
- const std::pair<typename CONTAINER::cell_iterator, unsigned int> pair2
- = std::make_pair(cell, 2*direction+1);
- pairs2.insert(pair2);
- }
- }
+ {
+ const typename CONTAINER::face_iterator face_1 = cell->face(2*direction);
+ const typename CONTAINER::face_iterator face_2 = cell->face(2*direction+1);
+
+ if (face_1->at_boundary() && face_1->boundary_indicator() == b_id)
+ {
+ const std::pair<typename CONTAINER::cell_iterator, unsigned int> pair1
+ = std::make_pair(cell, 2*direction);
+ pairs1.insert(pair1);
+ }
+
+ if (face_2->at_boundary() && face_2->boundary_indicator() == b_id)
+ {
+ const std::pair<typename CONTAINER::cell_iterator, unsigned int> pair2
+ = std::make_pair(cell, 2*direction+1);
+ pairs2.insert(pair2);
+ }
+ }
Assert (pairs1.size() == pairs2.size(),
ExcMessage ("Unmatched faces on periodic boundaries"));
#ifdef DEBUG
const unsigned int size_old = matched_pairs.size();
#endif
-
- // and call match_periodic_face_pairs that does the actual matching:
+
+ // and call match_periodic_face_pairs that does the actual matching:
match_periodic_face_pairs(pairs1, pairs2, direction, matched_pairs, offset);
#ifdef DEBUG
//check for standard orientation
const unsigned int size_new = matched_pairs.size();
for (unsigned int i = size_old; i < size_new; ++i)
- {
- Assert(matched_pairs[i].orientation == 1,
- ExcMessage("Found a face match with non standard orientation. "
- "This function is only suitable for meshes with cells "
- "in default orientation"));
- }
-#endif
+ {
+ Assert(matched_pairs[i].orientation == 1,
+ ExcMessage("Found a face match with non standard orientation. "
+ "This function is only suitable for meshes with cells "
+ "in default orientation"));
+ }
+#endif
}
if (points[n_intermediate_points].get() == 0)
{
std_cxx1x::shared_ptr<QGaussLobatto<1> >
- quadrature (new QGaussLobatto<1>(n_intermediate_points+2));
+ quadrature (new QGaussLobatto<1>(n_intermediate_points+2));
points[n_intermediate_points] = quadrature;
}
}
{
switch (dim)
{
- // in 2d, first check whether the two end points of the line are on the
- // axis of symmetry. if so, then return the mid point
+ // in 2d, first check whether the two end points of the line are on the
+ // axis of symmetry. if so, then return the mid point
case 2:
{
if ((line->vertex(0)(0) == this->center(0))
{
switch (dim)
{
- // in 2d, first check whether the two end points of the line are on the
- // axis of symmetry. if so, then return the mid point
+ // in 2d, first check whether the two end points of the line are on the
+ // axis of symmetry. if so, then return the mid point
case 2:
{
if ((line->vertex(0)(0) == this->center(0))
}
// first reserve, then resize. Otherwise the std library can decide to allocate
- // more entries.
- boundary_or_material_id.reserve (new_size);
- boundary_or_material_id.resize (new_size);
+ // more entries.
+ boundary_or_material_id.reserve (new_size);
+ boundary_or_material_id.resize (new_size);
user_data.reserve (new_size);
user_data.resize (new_size);
}
4*new_size-children.size(),
-1);
- // for the following two fields, we know exactly how many elements
- // we need, so first reserve then resize (resize itself, at least
- // with some compiler libraries, appears to round up the size it
- // actually reserves)
+ // for the following two fields, we know exactly how many elements
+ // we need, so first reserve then resize (resize itself, at least
+ // with some compiler libraries, appears to round up the size it
+ // actually reserves)
boundary_or_material_id.reserve (new_size);
boundary_or_material_id.resize (new_size);
active_fe_backup(dof_handler.levels.size ());
for (unsigned int level = 0; level<dof_handler.levels.size (); ++level)
std::swap (dof_handler.levels[level]->active_fe_indices,
- active_fe_backup[level]);
+ active_fe_backup[level]);
// delete all levels and set them up
// newly, since vectors are
// non-invalid value later on
dof_handler.faces->lines.dof_offsets
= std::vector<unsigned int> (dof_handler.tria->n_raw_lines(),
- (unsigned int)(-1));
+ (unsigned int)(-1));
dof_handler.faces->lines.dofs
= std::vector<types::global_dof_index> (n_line_slots,
DoFHandler<dim,spacedim>::invalid_dof_index);
{
dof_handler.faces->quads.dof_offsets
= std::vector<unsigned int>
- (dof_handler.tria->n_raw_quads(),
- (unsigned int)(-1));
+ (dof_handler.tria->n_raw_quads(),
+ (unsigned int)(-1));
dof_handler.faces->quads.dofs
= std::vector<types::global_dof_index> (n_quad_slots,
DoFHandler<dim,spacedim>::invalid_dof_index);
dofs_on_face.resize (dofs_per_face);
cell->face(f)->get_dof_indices (dofs_on_face,
- cell->active_fe_index());
+ cell->active_fe_index());
for (unsigned int i=0; i<dofs_per_face; ++i)
boundary_dofs.insert(dofs_on_face[i]);
};
// find the next cell active on this level
unsigned int next_cell = cell+1;
while ((next_cell<dof_offsets.size()) &&
- (dof_offsets[next_cell] == (offset_type)(-1)))
+ (dof_offsets[next_cell] == (offset_type)(-1)))
++next_cell;
const unsigned int next_offset = (next_cell < dof_offsets.size() ?
- dof_offsets[next_cell] :
- dof_indices.size());
+ dof_offsets[next_cell] :
+ dof_indices.size());
Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
ExcInternalError());
// find the next cell active on this level
unsigned int next_cell = cell+1;
while ((next_cell<dof_offsets.size()) &&
- (dof_offsets[next_cell] == (offset_type)(-1)))
+ (dof_offsets[next_cell] == (offset_type)(-1)))
++next_cell;
const unsigned int next_offset = (next_cell < dof_offsets.size() ?
- dof_offsets[next_cell] :
- dof_indices.size());
+ dof_offsets[next_cell] :
+ dof_indices.size());
Assert (next_offset-dof_offsets[cell] == fe_collection[active_fe_indices[cell]].template n_dofs_per_object<dim>(),
ExcInternalError());
// find the next cell active on this level
unsigned int next_cell = cell+1;
while ((next_cell<dof_offsets.size()) &&
- (dof_offsets[next_cell] == (offset_type)(-1)))
+ (dof_offsets[next_cell] == (offset_type)(-1)))
++next_cell;
const unsigned int next_offset = (next_cell < dof_offsets.size() ?
- dof_offsets[next_cell] :
- dof_indices.size());
+ dof_offsets[next_cell] :
+ dof_indices.size());
// set offset for this cell
new_dof_offsets[cell] = new_dof_indices.size();
DoFLevel::memory_consumption () const
{
return (MemoryConsumption::memory_consumption (active_fe_indices) +
- MemoryConsumption::memory_consumption (dof_indices) +
+ MemoryConsumption::memory_consumption (dof_indices) +
MemoryConsumption::memory_consumption (dof_offsets) +
MemoryConsumption::memory_consumption (cell_cache_offsets) +
MemoryConsumption::memory_consumption(cell_dof_indices_cache));
ConstraintMatrix::resolve_indices (std::vector<types::global_dof_index> &indices) const
{
const unsigned int indices_size = indices.size();
- const std::vector<std::pair<types::global_dof_index,double> > * line_ptr;
+ const std::vector<std::pair<types::global_dof_index,double> > *line_ptr;
for (unsigned int i=0; i<indices_size; ++i)
{
line_ptr = get_constraint_entries(indices[i]);
// transformation and target the wanted eigenvalues
if (transformation)
{
- // set transformation type if any
- transformation->set_context (solver_data->eps);
-
- // set target eigenvalues to solve for
- ierr = EPSSetTarget (solver_data->eps, target_eigenvalue);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+ // set transformation type if any
+ transformation->set_context (solver_data->eps);
+
+ // set target eigenvalues to solve for
+ ierr = EPSSetTarget (solver_data->eps, target_eigenvalue);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
}
// set which portion of the eigenspectrum to solve for
// get number of converged eigenstates
ierr = EPSGetConverged (solver_data->eps,
- reinterpret_cast<PetscInt*>(n_converged));
+ reinterpret_cast<PetscInt *>(n_converged));
AssertThrow (ierr == 0, ExcSLEPcError(ierr));
PetscInt n_iterations = 0;
// get the residual norm of the most extreme eigenvalue if and
// only if at least one eigenvector has converged.
if ((*n_converged)>0)
- {
- ierr = EPSComputeResidualNorm (solver_data->eps, 0, &residual_norm);
- AssertThrow (ierr == 0, ExcSLEPcError(ierr));
- }
+ {
+ ierr = EPSComputeResidualNorm (solver_data->eps, 0, &residual_norm);
+ AssertThrow (ierr == 0, ExcSLEPcError(ierr));
+ }
// check the solver state
const SolverControl::State state
const bool m_log_history,
const bool m_log_result)
:
- SolverControl (n, tolerance, m_log_history, m_log_result){}
+ SolverControl (n, tolerance, m_log_history, m_log_result) {}
IterationNumberControl::~IterationNumberControl()
template <class Matrix>
void SparseDirectMUMPS::initialize_matrix (const Matrix &matrix)
{
- Assert(matrix.n() == matrix.m(), ExcMessage("Matrix needs to be square."));
+ Assert(matrix.n() == matrix.m(), ExcMessage("Matrix needs to be square."));
// Check we haven't been here before:
}
-template <class CSP_t>
-void distribute_sparsity_pattern(CSP_t &csp,
+ template <class CSP_t>
+ void distribute_sparsity_pattern(CSP_t &csp,
const std::vector<IndexSet> &owned_set_per_cpu,
const MPI_Comm &mpi_comm,
const IndexSet &myrange)
-{
- size_type myid = Utilities::MPI::this_mpi_process(mpi_comm);
+ {
+ size_type myid = Utilities::MPI::this_mpi_process(mpi_comm);
- typedef std::map<size_type, std::vector<size_type> > map_vec_t;
- map_vec_t send_data;
+ typedef std::map<size_type, std::vector<size_type> > map_vec_t;
+ map_vec_t send_data;
- {
- unsigned int dest_cpu=0;
+ {
+ unsigned int dest_cpu=0;
- size_type n_local_rel_rows = myrange.n_elements();
- for (size_type row_idx=0; row_idx<n_local_rel_rows; ++row_idx)
- {
- size_type row=myrange.nth_index_in_set(row_idx);
+ size_type n_local_rel_rows = myrange.n_elements();
+ for (size_type row_idx=0; row_idx<n_local_rel_rows; ++row_idx)
+ {
+ size_type row=myrange.nth_index_in_set(row_idx);
- // calculate destination CPU, note that we start the search
- // at last destination cpu, because even if the owned ranges
- // are not contiguous, they hopefully consist of large blocks
- while (!owned_set_per_cpu[dest_cpu].is_element(row))
- {
- ++dest_cpu;
- if (dest_cpu==owned_set_per_cpu.size()) // wrap around
+ // calculate destination CPU, note that we start the search
+ // at last destination cpu, because even if the owned ranges
+ // are not contiguous, they hopefully consist of large blocks
+ while (!owned_set_per_cpu[dest_cpu].is_element(row))
+ {
+ ++dest_cpu;
+ if (dest_cpu==owned_set_per_cpu.size()) // wrap around
dest_cpu=0;
- }
+ }
- //skip myself
- if (dest_cpu==myid)
- continue;
+ //skip myself
+ if (dest_cpu==myid)
+ continue;
- size_type rlen = csp.row_length(row);
+ size_type rlen = csp.row_length(row);
- //skip empty lines
- if (!rlen)
- continue;
+ //skip empty lines
+ if (!rlen)
+ continue;
- //save entries
- std::vector<size_type> &dst = send_data[dest_cpu];
+ //save entries
+ std::vector<size_type> &dst = send_data[dest_cpu];
- dst.push_back(rlen); // number of entries
- dst.push_back(row); // row index
- for (size_type c=0; c<rlen; ++c)
- {
- //columns
- size_type column = csp.column_number(row, c);
- dst.push_back(column);
- }
- }
+ dst.push_back(rlen); // number of entries
+ dst.push_back(row); // row index
+ for (size_type c=0; c<rlen; ++c)
+ {
+ //columns
+ size_type column = csp.column_number(row, c);
+ dst.push_back(column);
+ }
+ }
- }
+ }
- unsigned int num_receive=0;
- {
- std::vector<unsigned int> send_to;
- send_to.reserve(send_data.size());
- for (map_vec_t::iterator it=send_data.begin(); it!=send_data.end(); ++it)
- send_to.push_back(it->first);
-
- num_receive =
- Utilities::MPI::
- compute_point_to_point_communication_pattern(mpi_comm, send_to).size();
- }
+ unsigned int num_receive=0;
+ {
+ std::vector<unsigned int> send_to;
+ send_to.reserve(send_data.size());
+ for (map_vec_t::iterator it=send_data.begin(); it!=send_data.end(); ++it)
+ send_to.push_back(it->first);
- std::vector<MPI_Request> requests(send_data.size());
+ num_receive =
+ Utilities::MPI::
+ compute_point_to_point_communication_pattern(mpi_comm, send_to).size();
+ }
+ std::vector<MPI_Request> requests(send_data.size());
- // send data
- {
- unsigned int idx=0;
- for (map_vec_t::iterator it=send_data.begin(); it!=send_data.end(); ++it, ++idx)
- MPI_Isend(&(it->second[0]),
- it->second.size(),
- DEAL_II_DOF_INDEX_MPI_TYPE,
- it->first,
- 124,
- mpi_comm,
- &requests[idx]);
- }
+
+ // send data
+ {
+ unsigned int idx=0;
+ for (map_vec_t::iterator it=send_data.begin(); it!=send_data.end(); ++it, ++idx)
+ MPI_Isend(&(it->second[0]),
+ it->second.size(),
+ DEAL_II_DOF_INDEX_MPI_TYPE,
+ it->first,
+ 124,
+ mpi_comm,
+ &requests[idx]);
+ }
//TODO: In the following, we read individual bytes and then reinterpret them
// as size_type objects. this is error prone. use properly typed reads that
// match the write above
- {
- //receive
- std::vector<size_type> recv_buf;
- for (unsigned int index=0; index<num_receive; ++index)
- {
- MPI_Status status;
- int len;
- MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, mpi_comm, &status);
- Assert (status.MPI_TAG==124, ExcInternalError());
+ {
+ //receive
+ std::vector<size_type> recv_buf;
+ for (unsigned int index=0; index<num_receive; ++index)
+ {
+ MPI_Status status;
+ int len;
+ MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, mpi_comm, &status);
+ Assert (status.MPI_TAG==124, ExcInternalError());
+
+ MPI_Get_count(&status, MPI_BYTE, &len);
+ Assert( len%sizeof(unsigned int)==0, ExcInternalError());
- MPI_Get_count(&status, MPI_BYTE, &len);
- Assert( len%sizeof(unsigned int)==0, ExcInternalError());
+ recv_buf.resize(len/sizeof(size_type));
- recv_buf.resize(len/sizeof(size_type));
+ MPI_Recv(&recv_buf[0], len, MPI_BYTE, status.MPI_SOURCE,
+ status.MPI_TAG, mpi_comm, &status);
- MPI_Recv(&recv_buf[0], len, MPI_BYTE, status.MPI_SOURCE,
- status.MPI_TAG, mpi_comm, &status);
+ std::vector<size_type>::const_iterator ptr = recv_buf.begin();
+ std::vector<size_type>::const_iterator end = recv_buf.end();
+ while (ptr+1<end)
+ {
+ size_type num=*(ptr++);
+ size_type row=*(ptr++);
+ for (unsigned int c=0; c<num; ++c)
+ {
+ csp.add(row, *ptr);
+ ptr++;
+ }
+ }
+ Assert(ptr==end, ExcInternalError());
+ }
+ }
- std::vector<size_type>::const_iterator ptr = recv_buf.begin();
- std::vector<size_type>::const_iterator end = recv_buf.end();
- while (ptr+1<end)
- {
- size_type num=*(ptr++);
- size_type row=*(ptr++);
- for (unsigned int c=0; c<num; ++c)
- {
- csp.add(row, *ptr);
- ptr++;
- }
- }
- Assert(ptr==end, ExcInternalError());
- }
+ // complete all sends, so that we can
+ // safely destroy the buffers.
+ MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
}
- // complete all sends, so that we can
- // safely destroy the buffers.
- MPI_Waitall(requests.size(), &requests[0], MPI_STATUSES_IGNORE);
-}
-
#endif
}
template void SparsityTools::distribute_sparsity_pattern
<BlockCompressedSimpleSparsityPattern>
(BlockCompressedSimpleSparsityPattern &csp,
- const std::vector<IndexSet> &owned_set_per_cpu,
- const MPI_Comm &mpi_comm,
- const IndexSet &myrange);
+ const std::vector<IndexSet> &owned_set_per_cpu,
+ const MPI_Comm &mpi_comm,
+ const IndexSet &myrange);
#endif
// we only allow assignment to vectors with the same number of blocks
// or to an empty BlockVector
Assert (n_blocks() == 0 || n_blocks() == v.n_blocks(),
- ExcDimensionMismatch(n_blocks(), v.n_blocks()));
+ ExcDimensionMismatch(n_blocks(), v.n_blocks()));
if (this->n_blocks() != v.n_blocks())
reinit(v.n_blocks());
/* -------------------------- PreconditionAMG -------------------------- */
- void
+ void
PreconditionIdentity::vmult(VectorBase &dst,
const VectorBase &src) const
{
dst = src;
}
- void
+ void
PreconditionIdentity::vmult(dealii::Vector<double> &dst,
const dealii::Vector<double> &src) const
{
Assert (false, ExcNotImplemented());
}
- // Introduce the preconditioner,
+ // Introduce the preconditioner,
// if the identity preconditioner is used,
// the precondioner is set to none, ...
if (preconditioner.preconditioner.use_count()!=0)
- {
- ierr = solver.SetPrecOperator (const_cast<Epetra_Operator *>
- (preconditioner.preconditioner.get()));
- AssertThrow (ierr == 0, ExcTrilinosError(ierr));
- }
+ {
+ ierr = solver.SetPrecOperator (const_cast<Epetra_Operator *>
+ (preconditioner.preconditioner.get()));
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
else
solver.SetAztecOption(AZ_precond,AZ_none);
:
(dim == 2 ?
(i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
+ :
+ (dim == 3 ?
+ (i<4*fe.dofs_per_vertex ?
+ i
:
- (dim == 3 ?
- (i<4*fe.dofs_per_vertex ?
- i
- :
- (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
- i+4*fe.dofs_per_vertex
- :
- i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+ (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
+ i+4*fe.dofs_per_vertex
:
- numbers::invalid_unsigned_int)));
+ i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
+ :
+ numbers::invalid_unsigned_int)));
Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
// make sure
const typename FunctionMap<dim>::type &function_map,
std::vector<IndexSet> &boundary_indices,
const ComponentMask &component_mask)
-{
+ {
Assert (boundary_indices.size() == dof.get_tria().n_global_levels(),
ExcDimensionMismatch (boundary_indices.size(),
dof.get_tria().n_global_levels()));
}
const unsigned int patch_idx =
- (*data.cell_to_patch_index_map)[cell_and_index->first->level()][cell_and_index->first->index()];
+ (*data.cell_to_patch_index_map)[cell_and_index->first->level()][cell_and_index->first->index()];
// did we mess up the indices?
Assert(patch_idx < patches.size(), ExcInternalError());
unsigned int active_index = 0;
cell_iterator cell = first_locally_owned_cell();
for (; cell != this->triangulation->end();
- cell = next_locally_owned_cell(cell))
+ cell = next_locally_owned_cell(cell))
{
// move forward until active_cell points at the cell (cell) we are looking
// at to compute the current active_index
while (active_cell!=this->triangulation->end()
- && cell->active()
- && active_cell_iterator(cell) != active_cell)
+ && cell->active()
+ && active_cell_iterator(cell) != active_cell)
{
++active_cell;
++active_index;
cell_to_patch_index_map[cell->level()].size(),
ExcInternalError());
Assert (active_index < this->triangulation->n_active_cells(),
- ExcInternalError());
+ ExcInternalError());
cell_to_patch_index_map[cell->level()][cell->index()] = all_cells.size();
all_cells.push_back (std::make_pair(cell, active_index));
std_cxx1x::bind(&DataOut<dim,DH>::build_one_patch,
this, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3,
curved_cell_region,std_cxx1x::ref(this->patches)),
- // no copy-local-to-global function needed here
- std_cxx1x::function<void (const ::dealii::DataOutBase::Patch<DH::dimension, DH::space_dimension> &)>(),
+ // no copy-local-to-global function needed here
+ std_cxx1x::function<void (const ::dealii::DataOutBase::Patch<DH::dimension, DH::space_dimension> &)>(),
thread_data,
sample_patch,
- // experimenting shows that we can make things run a bit
- // faster if we increase the number of cells we work on
- // per item (i.e., WorkStream's chunk_size argument,
- // about 10% improvement) and the items in flight at any
- // given time (another 5% on the testcase discussed in
- // @ref workstream_paper, on 32 cores) and if
- 8*multithread_info.n_threads(),
- 64);
+ // experimenting shows that we can make things run a bit
+ // faster if we increase the number of cells we work on
+ // per item (i.e., WorkStream's chunk_size argument,
+ // about 10% improvement) and the items in flight at any
+ // given time (another 5% on the testcase discussed in
+ // @ref workstream_paper, on 32 cores) and if
+ 8*multithread_info.n_threads(),
+ 64);
}
// convert the iterator to an active_iterator and advance this to the next
// active cell
typename Triangulation<DH::dimension,DH::space_dimension>::
- active_cell_iterator active_cell = cell;
+ active_cell_iterator active_cell = cell;
++active_cell;
return active_cell;
}
DataOut<dim,DH>::first_locally_owned_cell ()
{
typename DataOut<dim,DH>::cell_iterator
- cell = first_cell();
+ cell = first_cell();
// skip cells if the current one has no children (is active) and is a ghost
// or artificial cell
if (use_face_values == false)
{
dealii::hp::QCollection<dim>
- quadrature(QIterated<dim>(QTrapez<1>(), n_subdivisions));
+ quadrature(QIterated<dim>(QTrapez<1>(), n_subdivisions));
n_q_points = quadrature[0].size();
x_fe_values.resize(this->finite_elements.size());
for (unsigned int i=0; i<this->finite_elements.size(); ++i)
else
{
dealii::hp::QCollection<dim-1>
- quadrature(QIterated<dim-1>(QTrapez<1>(), n_subdivisions));
+ quadrature(QIterated<dim-1>(QTrapez<1>(), n_subdivisions));
n_q_points = quadrature[0].size();
x_fe_face_values.resize(this->finite_elements.size());
for (unsigned int i=0; i<this->finite_elements.size(); ++i)
{
Assert(data.x_fe_face_values.empty() == true, ExcInternalError());
dealii::hp::QCollection<dim>
- quadrature(QIterated<dim>(QTrapez<1>(), n_subdivisions));
+ quadrature(QIterated<dim>(QTrapez<1>(), n_subdivisions));
x_fe_values.resize(this->finite_elements.size());
for (unsigned int i=0; i<this->finite_elements.size(); ++i)
{
else
{
dealii::hp::QCollection<dim-1>
- quadrature(QIterated<dim-1>(QTrapez<1>(), n_subdivisions));
+ quadrature(QIterated<dim-1>(QTrapez<1>(), n_subdivisions));
x_fe_face_values.resize(this->finite_elements.size());
for (unsigned int i=0; i<this->finite_elements.size(); ++i)
{
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
DataOut_DoFData<DH,patch_dim,patch_space_dim>::DataOut_DoFData ()
:
triangulation(0,typeid(*this).name()),
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
template <class VECTOR>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
template <class VECTOR>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
template <class VECTOR>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
template <class VECTOR>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
template <class VECTOR>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
template <class VECTOR>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
void DataOut_DoFData<DH,patch_dim,patch_space_dim>::clear_data_vectors ()
{
dof_data.erase (dof_data.begin(), dof_data.end());
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
clear_input_data_references ()
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
void
DataOut_DoFData<DH,patch_dim,patch_space_dim>::clear ()
{
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
std::vector<std::string>
DataOut_DoFData<DH,patch_dim,patch_space_dim>::
get_dataset_names () const
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
std::vector<std_cxx1x::tuple<unsigned int, unsigned int, std::string> >
DataOut_DoFData<DH,patch_dim,patch_space_dim>::get_vector_data_ranges () const
{
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
const std::vector< dealii::DataOutBase::Patch<patch_dim, patch_space_dim> > &
DataOut_DoFData<DH,patch_dim,patch_space_dim>::get_patches () const
{
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
std::vector<std_cxx1x::shared_ptr<dealii::hp::FECollection<DH::dimension,DH::space_dimension> > >
DataOut_DoFData<DH,patch_dim,patch_space_dim>::get_finite_elements() const
{
const unsigned int dhdim = DH::dimension;
const unsigned int dhspacedim = DH::space_dimension;
std::vector<std_cxx1x::shared_ptr<dealii::hp::FECollection<dhdim,dhspacedim> > >
- finite_elements(this->dof_data.size());
+ finite_elements(this->dof_data.size());
for (unsigned int i=0; i<this->dof_data.size(); ++i)
{
Assert (dof_data[i]->dof_handler != 0, ExcNoDoFHandlerSelected ());
template <class DH,
- int patch_dim, int patch_space_dim>
+ int patch_dim, int patch_space_dim>
std::size_t
DataOut_DoFData<DH,patch_dim,patch_space_dim>::memory_consumption () const
{
// first fill dof_data
for (unsigned int dataset=0; dataset<this->dof_data.size(); ++dataset)
{
- const FEValuesBase<DH::dimension> & fe_patch_values
+ const FEValuesBase<DH::dimension> &fe_patch_values
= data.get_present_fe_values(dataset);
const unsigned int n_components
= fe_patch_values.get_fe().n_components();
template <class DerivativeDescription, int dim,
- template <int, int> class DH, class InputVector, int spacedim>
+ template <int, int> class DH, class InputVector, int spacedim>
void
DerivativeApproximation::
approximate_derivative (const Mapping<dim,spacedim> &mapping,
// Only act on the locally owned cells
typedef FilteredIterator<typename DH<dim,spacedim>::active_cell_iterator> CellFilter;
-
+
typedef std_cxx1x::tuple<CellFilter,Vector<float>::iterator>
- Iterators;
+ Iterators;
SynchronousIterators<Iterators> begin(Iterators (CellFilter(IteratorFilters::LocallyOwnedCell(),
- dof_handler.begin_active()),derivative_norm.begin())),
- end(Iterators (CellFilter(IteratorFilters::LocallyOwnedCell(),dof_handler.end()),
- derivative_norm.end()));
-
- // There is no need for a copier because there is no conflict between threads
+ dof_handler.begin_active()),derivative_norm.begin())),
+ end(Iterators (CellFilter(IteratorFilters::LocallyOwnedCell(),dof_handler.end()),
+ derivative_norm.end()));
+
+ // There is no need for a copier because there is no conflict between threads
// to write in derivative_norm. Scratch and CopyData are also useless.
WorkStream::run(begin,end,
- static_cast<std_cxx1x::function<void (SynchronousIterators<Iterators> const&,
- internal::Assembler::Scratch const&,internal::Assembler::CopyData &)> >
- (std_cxx1x::bind(&DerivativeApproximation::template approximate<DerivativeDescription,dim,DH,
- InputVector,spacedim>,
- std_cxx1x::_1,
- std_cxx1x::cref(mapping),
- std_cxx1x::cref(dof_handler),
- std_cxx1x::cref(solution),component)),
- std_cxx1x::function<void (internal::Assembler::CopyData const &)> (),
- internal::Assembler::Scratch (),internal::Assembler::CopyData ());
+ static_cast<std_cxx1x::function<void (SynchronousIterators<Iterators> const &,
+ internal::Assembler::Scratch const &,internal::Assembler::CopyData &)> >
+ (std_cxx1x::bind(&DerivativeApproximation::template approximate<DerivativeDescription,dim,DH,
+ InputVector,spacedim>,
+ std_cxx1x::_1,
+ std_cxx1x::cref(mapping),
+ std_cxx1x::cref(dof_handler),
+ std_cxx1x::cref(solution),component)),
+ std_cxx1x::function<void (internal::Assembler::CopyData const &)> (),
+ internal::Assembler::Scratch (),internal::Assembler::CopyData ());
}
template <class DerivativeDescription, int dim,
- template <int, int> class DH, class InputVector, int spacedim>
+ template <int, int> class DH, class InputVector, int spacedim>
void
DerivativeApproximation::approximate (SynchronousIterators<std_cxx1x::tuple<FilteredIterator<typename DH<dim,spacedim>::active_cell_iterator>,Vector<float>::iterator> > const &cell,
const Mapping<dim,spacedim> &mapping,
// call the function doing the actual
// work on this cell
DerivativeApproximation::template approximate_cell<DerivativeDescription,dim,DH,InputVector>
- (mapping,dof_handler,solution,component,std_cxx1x::get<0>(cell.iterators),derivative);
+ (mapping,dof_handler,solution,component,std_cxx1x::get<0>(cell.iterators),derivative);
// evaluate the norm and fill the vector
//*derivative_norm_on_this_cell
*std_cxx1x::get<1>(cell.iterators) = DerivativeDescription::derivative_norm (derivative);
template <class DerivativeDescription, int dim,
- template <int, int> class DH, class InputVector, int spacedim>
+ template <int, int> class DH, class InputVector, int spacedim>
void
DerivativeApproximation::
approximate_cell (const Mapping<dim,spacedim> &mapping,
namespace AssemblerData
{
template <int dim,
- int spacedim>
+ int spacedim>
struct Scratch
{
Scratch (const ::dealii::hp::FECollection<dim,spacedim> &fe,
template <int dim,
- int spacedim,
- typename CellIterator>
+ int spacedim,
+ typename CellIterator>
void mass_assembler (const CellIterator &cell,
MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data,
MatrixCreator::internal::AssemblerData::CopyData ©_data)
template <int dim,
- int spacedim,
- typename CellIterator>
+ int spacedim,
+ typename CellIterator>
void laplace_assembler (const CellIterator &cell,
MatrixCreator::internal::AssemblerData::Scratch<dim,spacedim> &data,
MatrixCreator::internal::AssemblerData::CopyData ©_data)
template <typename MatrixType,
- typename VectorType>
+ typename VectorType>
void copy_local_to_global (const AssemblerData::CopyData &data,
MatrixType *matrix,
VectorType *right_hand_side)
template <typename DH>
CopyData<DH>::CopyData(CopyData const &data) :
- dofs_per_cell(data.dofs_per_cell),
- dofs(data.dofs),
- dof_is_on_face(data.dof_is_on_face),
- cell(data.cell),
- cell_matrix(data.cell_matrix),
- cell_vector(data.cell_vector)
+ dofs_per_cell(data.dofs_per_cell),
+ dofs(data.dofs),
+ dof_is_on_face(data.dof_is_on_face),
+ cell(data.cell),
+ cell_matrix(data.cell_matrix),
+ cell_vector(data.cell_vector)
{}
}
}
typename FunctionMap<spacedim>::type const &boundary_functions,
Function<spacedim> const *const coefficient,
std::vector<unsigned int> const &component_mapping)
-
+
{
// All assertions for this function
// are in the calling function
const bool fe_is_primitive = fe.is_primitive();
const unsigned int dofs_per_face = fe.dofs_per_face;
-
+
copy_data.cell = cell;
copy_data.dofs_per_cell = fe.dofs_per_cell;
std::vector<double> rhs_values_scalar (fe_values.n_quadrature_points);
std::vector<Vector<double> > rhs_values_system (fe_values.n_quadrature_points,
- Vector<double>(n_function_components));
+ Vector<double>(n_function_components));
copy_data.dofs.resize(copy_data.dofs_per_cell);
cell->get_dof_indices (copy_data.dofs);
boundary_functions.end())
{
copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
- copy_data.dofs_per_cell));
+ copy_data.dofs_per_cell));
copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
fe_values.reinit (cell, face);
}
}
copy_data.cell_vector.back()(i) += fe_values.shape_value(i,point)
- * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first])
- * weight;
+ * rhs_values_system[point](component_mapping[fe.system_to_component_index(i).first])
+ * weight;
}
else
{
* normal_adjustment[point][comp]
* weight * coefficient_vector_values[point](comp);
copy_data.cell_vector.back()(i) += fe_values.shape_value_component(i,point,comp) *
- rhs_values_system[point](component_mapping[comp])
- * normal_adjustment[point][comp]
- * weight;
+ rhs_values_system[point](component_mapping[comp])
+ * normal_adjustment[point][comp]
+ * weight;
}
}
}
copy_data.dofs[i])
!=
dofs_on_face_vector.end());
- }
+ }
}
template <int dim,int spacedim>
// searches.
unsigned int pos(0);
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- {
- // check if this face is on that part of
- // the boundary we are interested in
- if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) !=
- boundary_functions.end())
- {
- for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
- {
- if (copy_data.dof_is_on_face[pos][i] &&
- dof_to_boundary_mapping[copy_data.dofs[i]] != numbers::invalid_dof_index)
- {
- for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
- if (copy_data.dof_is_on_face[pos][j] &&
- dof_to_boundary_mapping[copy_data.dofs[j]] != numbers::invalid_dof_index)
- {
- Assert(numbers::is_finite(copy_data.cell_matrix[pos](i,j)),
- ExcNumberNotFinite());
- matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
- dof_to_boundary_mapping[copy_data.dofs[j]],
- copy_data.cell_matrix[pos](i,j));
- }
- Assert(numbers::is_finite(copy_data.cell_vector[pos](i)), ExcNumberNotFinite());
- rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += copy_data.cell_vector[pos](i);
- }
- }
- ++pos;
- }
- }
+ {
+ // check if this face is on that part of
+ // the boundary we are interested in
+ if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) !=
+ boundary_functions.end())
+ {
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ {
+ if (copy_data.dof_is_on_face[pos][i] &&
+ dof_to_boundary_mapping[copy_data.dofs[i]] != numbers::invalid_dof_index)
+ {
+ for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+ if (copy_data.dof_is_on_face[pos][j] &&
+ dof_to_boundary_mapping[copy_data.dofs[j]] != numbers::invalid_dof_index)
+ {
+ Assert(numbers::is_finite(copy_data.cell_matrix[pos](i,j)),
+ ExcNumberNotFinite());
+ matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
+ dof_to_boundary_mapping[copy_data.dofs[j]],
+ copy_data.cell_matrix[pos](i,j));
+ }
+ Assert(numbers::is_finite(copy_data.cell_vector[pos](i)), ExcNumberNotFinite());
+ rhs_vector(dof_to_boundary_mapping[copy_data.dofs[i]]) += copy_data.cell_vector[pos](i);
+ }
+ }
+ ++pos;
+ }
+ }
}
template <>
void
create_boundary_mass_matrix_1<2,3> (DoFHandler<2,3>::active_cell_iterator const &cell,
- MatrixCreator::internal::AssemblerBoundary::Scratch const
- &scratch,
+ MatrixCreator::internal::AssemblerBoundary::Scratch const
+ &scratch,
MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<2,
- 3> > ©_data,
- Mapping<2,3> const &mapping,
- FiniteElement<2,3> const &fe,
- Quadrature<1> const &q,
- FunctionMap<3>::type const &boundary_functions,
- Function<3> const *const coefficient,
- std::vector<unsigned int> const &component_mapping)
+ 3> > ©_data,
+ Mapping<2,3> const &mapping,
+ FiniteElement<2,3> const &fe,
+ Quadrature<1> const &q,
+ FunctionMap<3>::type const &boundary_functions,
+ Function<3> const *const coefficient,
+ std::vector<unsigned int> const &component_mapping)
{
Assert(false,ExcNotImplemented());
}
template <>
void
create_boundary_mass_matrix_1<1,3> (DoFHandler<1,3>::active_cell_iterator const &cell,
- MatrixCreator::internal::AssemblerBoundary::Scratch const
+ MatrixCreator::internal::AssemblerBoundary::Scratch const
&scratch,
MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<1,
- 3> > ©_data,
+ 3> > ©_data,
Mapping<1,3> const &mapping,
FiniteElement<1,3> const &fe,
Quadrature<0> const &q,
FunctionMap<3>::type const &boundary_functions,
Function<3> const *const coefficient,
- std::vector<unsigned int> const &component_mapping)
+ std::vector<unsigned int> const &component_mapping)
{
Assert(false,ExcNotImplemented());
}
MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > copy_data;
WorkStream::run(dof.begin_active(),dof.end(),
- static_cast<std_cxx1x::function<void (typename DoFHandler<dim,spacedim>::active_cell_iterator
- const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
- MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > &)> >
- (std_cxx1x::bind(&create_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,std_cxx1x::_2,
- std_cxx1x::_3,
- std_cxx1x::cref(mapping),std_cxx1x::cref(fe),std_cxx1x::cref(q),
- std_cxx1x::cref(boundary_functions),coefficient,
- std_cxx1x::cref(component_mapping))),
- static_cast<std_cxx1x::function<void (MatrixCreator::internal::AssemblerBoundary
- ::CopyData<DoFHandler<dim,spacedim> > const &)> > (std_cxx1x::bind(
- ©_boundary_mass_matrix_1<dim,spacedim>,
- std_cxx1x::_1,
- std_cxx1x::cref(boundary_functions),
- std_cxx1x::cref(dof_to_boundary_mapping),
- std_cxx1x::ref(matrix),
- std_cxx1x::ref(rhs_vector))),
- scratch,
- copy_data);
+ static_cast<std_cxx1x::function<void (typename DoFHandler<dim,spacedim>::active_cell_iterator
+ const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<DoFHandler<dim,spacedim> > &)> >
+ (std_cxx1x::bind(&create_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,std_cxx1x::_2,
+ std_cxx1x::_3,
+ std_cxx1x::cref(mapping),std_cxx1x::cref(fe),std_cxx1x::cref(q),
+ std_cxx1x::cref(boundary_functions),coefficient,
+ std_cxx1x::cref(component_mapping))),
+ static_cast<std_cxx1x::function<void (MatrixCreator::internal::AssemblerBoundary
+ ::CopyData<DoFHandler<dim,spacedim> > const &)> > (std_cxx1x::bind(
+ ©_boundary_mass_matrix_1<dim,spacedim>,
+ std_cxx1x::_1,
+ std_cxx1x::cref(boundary_functions),
+ std_cxx1x::cref(dof_to_boundary_mapping),
+ std_cxx1x::ref(matrix),
+ std_cxx1x::ref(rhs_vector))),
+ scratch,
+ copy_data);
}
template <int dim, int spacedim>
void
- create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator const
+ create_hp_boundary_mass_matrix_1 (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator const
&cell,
MatrixCreator::internal::AssemblerBoundary::Scratch const &scratch,
MatrixCreator::internal::AssemblerBoundary
const FEFaceValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values ();
copy_data.cell_matrix.push_back(FullMatrix<double> (copy_data.dofs_per_cell,
- copy_data.dofs_per_cell));
+ copy_data.dofs_per_cell));
copy_data.cell_vector.push_back(Vector<double> (copy_data.dofs_per_cell));
if (fe_is_system)
}
copy_data.cell_vector.back()(i) += v *
- rhs_values_system[point](
- component_mapping[fe.system_to_component_index(i).first]) * weight;
+ rhs_values_system[point](
+ component_mapping[fe.system_to_component_index(i).first]) * weight;
}
}
}
copy_data.cell_matrix.back()(i,j) +=
(u * v * weight * coefficient_vector_values[point](component_i));
}
- copy_data.cell_vector.back()(i) += v *
- rhs_values_system[point](component_mapping[component_i]) * weight;
+ copy_data.cell_vector.back()(i) += v *
+ rhs_values_system[point](component_mapping[component_i]) * weight;
}
}
}
copy_data.cell_matrix.back()(i,j) += (u * v * weight);
}
copy_data.cell_vector.back()(i) += v *
- rhs_values_system[point](
- fe.system_to_component_index(i).first) *
- weight;
+ rhs_values_system[point](
+ fe.system_to_component_index(i).first) *
+ weight;
}
}
}
for (unsigned int j=0; j<fe_values.dofs_per_cell; ++j)
{
const double u = fe_values.shape_value(j,point);
- copy_data.cell_matrix.back()(i,j) += (u * v * weight *
- coefficient_values[point]);
+ copy_data.cell_matrix.back()(i,j) += (u * v * weight *
+ coefficient_values[point]);
}
copy_data.cell_vector.back()(i) += v * rhs_values_scalar[point] *weight;
}
// searches.
unsigned int pos(0);
for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
- {
- // check if this face is on that part of
- // the boundary we are interested in
- if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) !=
- boundary_functions.end())
- {
+ {
+ // check if this face is on that part of
+ // the boundary we are interested in
+ if (boundary_functions.find(copy_data.cell->face(face)->boundary_indicator()) !=
+ boundary_functions.end())
+ {
#ifdef DEBUG
- // in debug mode: compute an element in the matrix which is
- // guaranteed to belong to a boundary dof. We do this to check that the
- // entries in the cell matrix are guaranteed to be zero if the
- // respective dof is not on the boundary. Since because of
- // round-off, the actual value of the matrix entry may be
- // only close to zero, we assert that it is small relative to an element
- // which is guaranteed to be nonzero. (absolute smallness does not
- // suffice since the size of the domain scales in here)
- //
- // for this purpose we seek the diagonal of the matrix, where there
- // must be an element belonging to the boundary. we take the maximum
- // diagonal entry.
- types::global_dof_index max_element = static_cast<types::global_dof_index>(0);
- for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
- i!=dof_to_boundary_mapping.end(); ++i)
- if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) &&
- (*i > max_element))
- max_element = *i;
- Assert (max_element == matrix.n()-1, ExcInternalError());
-
- double max_diag_entry = 0;
- for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
- if (std::fabs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry)
- max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i));
+ // in debug mode: compute an element in the matrix which is
+ // guaranteed to belong to a boundary dof. We do this to check that the
+ // entries in the cell matrix are guaranteed to be zero if the
+ // respective dof is not on the boundary. Since because of
+ // round-off, the actual value of the matrix entry may be
+ // only close to zero, we assert that it is small relative to an element
+ // which is guaranteed to be nonzero. (absolute smallness does not
+ // suffice since the size of the domain scales in here)
+ //
+ // for this purpose we seek the diagonal of the matrix, where there
+ // must be an element belonging to the boundary. we take the maximum
+ // diagonal entry.
+ types::global_dof_index max_element = static_cast<types::global_dof_index>(0);
+ for (std::vector<types::global_dof_index>::const_iterator i=dof_to_boundary_mapping.begin();
+ i!=dof_to_boundary_mapping.end(); ++i)
+ if ((*i != hp::DoFHandler<dim,spacedim>::invalid_dof_index) &&
+ (*i > max_element))
+ max_element = *i;
+ Assert (max_element == matrix.n()-1, ExcInternalError());
+
+ double max_diag_entry = 0;
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ if (std::fabs(copy_data.cell_matrix[pos](i,i)) > max_diag_entry)
+ max_diag_entry = std::fabs(copy_data.cell_matrix[pos](i,i));
#endif
- for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ for (unsigned int i=0; i<copy_data.dofs_per_cell; ++i)
+ for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
+ {
+ if (copy_data.dof_is_on_face[pos][i] && copy_data.dof_is_on_face[pos][j])
+ matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
+ dof_to_boundary_mapping[copy_data.dofs[j]],
+ copy_data.cell_matrix[pos](i,j));
+ else
+ {
+ // assume that all shape functions that are nonzero on the boundary
+ // are also listed in the @p{dof_to_boundary} mapping. if that
+ // is not the case, then the boundary mass matrix does not
+ // make that much sense anyway, as it only contains entries for
+ // parts of the functions living on the boundary
+ //
+ // these, we may compare here for relative smallness of all
+ // entries in the local matrix which are not taken over to
+ // the global one
+ Assert (std::fabs(copy_data.cell_matrix[pos](i,j)) <= 1e-10 * max_diag_entry,
+ ExcInternalError ());
+ }
+ }
+
for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
- {
- if (copy_data.dof_is_on_face[pos][i] && copy_data.dof_is_on_face[pos][j])
- matrix.add(dof_to_boundary_mapping[copy_data.dofs[i]],
- dof_to_boundary_mapping[copy_data.dofs[j]],
- copy_data.cell_matrix[pos](i,j));
+ if (copy_data.dof_is_on_face[pos][j])
+ rhs_vector(dof_to_boundary_mapping[copy_data.dofs[j]]) += copy_data.cell_vector[pos](j);
else
{
- // assume that all shape functions that are nonzero on the boundary
- // are also listed in the @p{dof_to_boundary} mapping. if that
- // is not the case, then the boundary mass matrix does not
- // make that much sense anyway, as it only contains entries for
- // parts of the functions living on the boundary
- //
- // these, we may compare here for relative smallness of all
- // entries in the local matrix which are not taken over to
- // the global one
- Assert (std::fabs(copy_data.cell_matrix[pos](i,j)) <= 1e-10 * max_diag_entry,
- ExcInternalError ());
+ // compare here for relative
+ // smallness
+ Assert (std::fabs(copy_data.cell_vector[pos](j)) <= 1e-10 * max_diag_entry,
+ ExcInternalError());
}
- }
-
- for (unsigned int j=0; j<copy_data.dofs_per_cell; ++j)
- if (copy_data.dof_is_on_face[pos][j])
- rhs_vector(dof_to_boundary_mapping[copy_data.dofs[j]]) += copy_data.cell_vector[pos](j);
- else
- {
- // compare here for relative
- // smallness
- Assert (std::fabs(copy_data.cell_vector[pos](j)) <= 1e-10 * max_diag_entry,
- ExcInternalError());
- }
- ++pos;
- }
- }
+ ++pos;
+ }
+ }
}
}
MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > copy_data;
WorkStream::run(dof.begin_active(),dof.end(),
- static_cast<std_cxx1x::function<void (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
- const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
- MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > &)> >
- (std_cxx1x::bind( &create_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,std_cxx1x::_2,
- std_cxx1x::_3,
- std_cxx1x::cref(mapping),std_cxx1x::cref(fe_collection),std_cxx1x::cref(q),
- std_cxx1x::cref(boundary_functions),coefficient,
- std_cxx1x::cref(component_mapping))),
- static_cast<std_cxx1x::function<void (MatrixCreator::internal::AssemblerBoundary
- ::CopyData<hp::DoFHandler<dim,spacedim> > const &)> > (
- std_cxx1x::bind( ©_hp_boundary_mass_matrix_1<dim,spacedim>,
- std_cxx1x::_1,
- std_cxx1x::cref(boundary_functions),
- std_cxx1x::cref(dof_to_boundary_mapping),
- std_cxx1x::ref(matrix),
- std_cxx1x::ref(rhs_vector))),
- scratch,
- copy_data);
+ static_cast<std_cxx1x::function<void (typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
+ const &,MatrixCreator::internal::AssemblerBoundary::Scratch const &,
+ MatrixCreator::internal::AssemblerBoundary::CopyData<hp::DoFHandler<dim,spacedim> > &)> >
+ (std_cxx1x::bind( &create_hp_boundary_mass_matrix_1<dim,spacedim>,std_cxx1x::_1,std_cxx1x::_2,
+ std_cxx1x::_3,
+ std_cxx1x::cref(mapping),std_cxx1x::cref(fe_collection),std_cxx1x::cref(q),
+ std_cxx1x::cref(boundary_functions),coefficient,
+ std_cxx1x::cref(component_mapping))),
+ static_cast<std_cxx1x::function<void (MatrixCreator::internal::AssemblerBoundary
+ ::CopyData<hp::DoFHandler<dim,spacedim> > const &)> > (
+ std_cxx1x::bind( ©_hp_boundary_mass_matrix_1<dim,spacedim>,
+ std_cxx1x::_1,
+ std_cxx1x::cref(boundary_functions),
+ std_cxx1x::cref(dof_to_boundary_mapping),
+ std_cxx1x::ref(matrix),
+ std_cxx1x::ref(rhs_vector))),
+ scratch,
+ copy_data);
}
static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
&MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>,
std_cxx1x::bind (&MatrixCreator::internal::
- copy_local_to_global<SparseMatrix<double>, Vector<double> >,
- std_cxx1x::_1,
- &matrix,
- (Vector<double> *)NULL),
+ copy_local_to_global<SparseMatrix<double>, Vector<double> >,
+ std_cxx1x::_1,
+ &matrix,
+ (Vector<double> *)NULL),
assembler_data,
copy_data);
}
static_cast<typename DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
&MatrixCreator::internal::laplace_assembler<dim, spacedim, typename DoFHandler<dim,spacedim>::active_cell_iterator>,
std_cxx1x::bind (&MatrixCreator::internal::
- copy_local_to_global<SparseMatrix<double>, Vector<double> >,
- std_cxx1x::_1,
- &matrix,
- &rhs_vector),
+ copy_local_to_global<SparseMatrix<double>, Vector<double> >,
+ std_cxx1x::_1,
+ &matrix,
+ &rhs_vector),
assembler_data,
copy_data);
}
static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
&MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>,
std_cxx1x::bind (&MatrixCreator::internal::
- copy_local_to_global<SparseMatrix<double>, Vector<double> >,
- std_cxx1x::_1,
- &matrix,
- (Vector<double> *)0),
+ copy_local_to_global<SparseMatrix<double>, Vector<double> >,
+ std_cxx1x::_1,
+ &matrix,
+ (Vector<double> *)0),
assembler_data,
copy_data);
}
static_cast<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>(dof.end()),
&MatrixCreator::internal::laplace_assembler<dim, spacedim, typename hp::DoFHandler<dim,spacedim>::active_cell_iterator>,
std_cxx1x::bind (&MatrixCreator::internal::
- copy_local_to_global<SparseMatrix<double>, Vector<double> >,
- std_cxx1x::_1,
- &matrix,
- &rhs_vector),
+ copy_local_to_global<SparseMatrix<double>, Vector<double> >,
+ std_cxx1x::_1,
+ &matrix,
+ &rhs_vector),
assembler_data,
copy_data);
}
local_values.reinit (cell->has_children() ?
cell->child(0)->get_fe().dofs_per_cell
: cell->get_fe().dofs_per_cell, true);
- // do the interpolation. we get into trouble if the
- // interpolation_hp(new,old) matrix hasn't been computed.
- // this can happen if the respective elements don't support
- // the corresponding interpolation; if that's the case, then
- // the computation of the matrix simply sets the matrix
- // back to size zero. so if we get here and that is
- // the wrong size, then this may be because the elements
- // haven't implemented the correct function yet
- //
- // there is one wrinkle. we would like to only error out if
- // the size of the matrix is 0 times 0 but at least one
- // of the elements has more than one dof per cell. the
- // problem is that if you reinit a matrix to 4x0, it automatically
- // sets its size to 0x0. so we can only execute the following
- // test if *both* elements have dofs_per_cell>0, not if *at
- // least one* have.
- Assert (! ((interpolation_hp(new_fe_index,old_index).m() == 0)
- &&
- (interpolation_hp(new_fe_index,old_index).n() == 0)
- &&
- ((dof_handler->get_fe()[new_fe_index].dofs_per_cell > 0)
- &&
- (dof_handler->get_fe()[old_index].dofs_per_cell > 0))),
- ExcMessage ("The interpolation between two different "
- "elements you are trying to use here has "
- "not been implemented for this pair of "
- "elements!"));
+ // do the interpolation. we get into trouble if the
+ // interpolation_hp(new,old) matrix hasn't been computed.
+ // this can happen if the respective elements don't support
+ // the corresponding interpolation; if that's the case, then
+ // the computation of the matrix simply sets the matrix
+ // back to size zero. so if we get here and that is
+ // the wrong size, then this may be because the elements
+ // haven't implemented the correct function yet
+ //
+ // there is one wrinkle. we would like to only error out if
+ // the size of the matrix is 0 times 0 but at least one
+ // of the elements has more than one dof per cell. the
+ // problem is that if you reinit a matrix to 4x0, it automatically
+ // sets its size to 0x0. so we can only execute the following
+ // test if *both* elements have dofs_per_cell>0, not if *at
+ // least one* have.
+ Assert (! ((interpolation_hp(new_fe_index,old_index).m() == 0)
+ &&
+ (interpolation_hp(new_fe_index,old_index).n() == 0)
+ &&
+ ((dof_handler->get_fe()[new_fe_index].dofs_per_cell > 0)
+ &&
+ (dof_handler->get_fe()[old_index].dofs_per_cell > 0))),
+ ExcMessage ("The interpolation between two different "
+ "elements you are trying to use here has "
+ "not been implemented for this pair of "
+ "elements!"));
// simple case where all children have the
// same FE index: just interpolate to their FE
// first and then use the standard routines
- if (tmp.size() > 0)
- interpolation_hp(new_fe_index,old_index).vmult (local_values, tmp);
- else
- local_values = 0;
+ if (tmp.size() > 0)
+ interpolation_hp(new_fe_index,old_index).vmult (local_values, tmp);
+ else
+ local_values = 0;
}
if (cell->has_children() == false)
interpolation_hp(c_index,old_index).n());
local_values.reinit(cell->child(child)->get_fe().dofs_per_cell, true);
- // do the interpolation. same problem as above
- Assert (! ((interpolation_hp(c_index,old_index).m() == 0)
- &&
- (interpolation_hp(c_index,old_index).n() == 0)
- &&
- ((dof_handler->get_fe()[c_index].dofs_per_cell > 0)
- &&
- (dof_handler->get_fe()[old_index].dofs_per_cell > 0))),
- ExcMessage ("The interpolation between two different "
- "elements you are trying to use here has "
- "not been implemented for this pair of "
- "elements!"));
-
- if (tmp.size() > 0)
- interpolation_hp(c_index,old_index).vmult (local_values, tmp);
+ // do the interpolation. same problem as above
+ Assert (! ((interpolation_hp(c_index,old_index).m() == 0)
+ &&
+ (interpolation_hp(c_index,old_index).n() == 0)
+ &&
+ ((dof_handler->get_fe()[c_index].dofs_per_cell > 0)
+ &&
+ (dof_handler->get_fe()[old_index].dofs_per_cell > 0))),
+ ExcMessage ("The interpolation between two different "
+ "elements you are trying to use here has "
+ "not been implemented for this pair of "
+ "elements!"));
+
+ if (tmp.size() > 0)
+ interpolation_hp(c_index,old_index).vmult (local_values, tmp);
else
local_values = 0;
}
projects / dealii-svn.git / commitdiff