</p>
<ol>
- <li> Changed: ConstraintMatrix::distribute_local_to_global() now requires
- matching data types. This is done to correctly handle complex-valued case.
+ <li> Changed: ConstraintMatrix::distribute_local_to_global() and numerous
+ functions in VectorTools namespace now require matching data types.
+ This is done to correctly handle complex-valued case.
<br>
(Denis Davydov, 2016/02/22)
</li>
template <typename VectorType, int dim, int spacedim, template <int, int> class DoFHandlerType>
void interpolate (const Mapping<dim,spacedim> &mapping,
const DoFHandlerType<dim,spacedim> &dof,
- const Function<spacedim,double> &function,
+ const Function<spacedim,typename VectorType::value_type> &function,
VectorType &vec);
/**
*/
template <typename VectorType, typename DoFHandlerType>
void interpolate (const DoFHandlerType &dof,
- const Function<DoFHandlerType::space_dimension,double> &function,
+ const Function<DoFHandlerType::space_dimension,typename VectorType::value_type> &function,
VectorType &vec);
/**
interpolate_based_on_material_id
(const Mapping<DoFHandlerType::dimension, DoFHandlerType::space_dimension> &mapping,
const DoFHandlerType &dof_handler,
- const std::map<types::material_id, const Function<DoFHandlerType::space_dimension, double> *> &function_map,
+ const std::map<types::material_id, const Function<DoFHandlerType::space_dimension, typename VectorType::value_type> *> &function_map,
VectorType &dst,
const ComponentMask &component_mask = ComponentMask());
const DoFHandler<dim,spacedim> &dof,
const ConstraintMatrix &constraints,
const Quadrature<dim> &quadrature,
- const Function<spacedim,double> &function,
+ const Function<spacedim,typename VectorType::value_type> &function,
VectorType &vec,
const bool enforce_zero_boundary = false,
const Quadrature<dim-1> &q_boundary = (dim > 1 ?
void project (const DoFHandler<dim,spacedim> &dof,
const ConstraintMatrix &constraints,
const Quadrature<dim> &quadrature,
- const Function<spacedim,double> &function,
+ const Function<spacedim,typename VectorType::value_type> &function,
VectorType &vec,
const bool enforce_zero_boundary = false,
const Quadrature<dim-1> &q_boundary = (dim > 1 ?
const hp::DoFHandler<dim,spacedim> &dof,
const ConstraintMatrix &constraints,
const hp::QCollection<dim> &quadrature,
- const Function<spacedim,double> &function,
+ const Function<spacedim,typename VectorType::value_type> &function,
VectorType &vec,
const bool enforce_zero_boundary = false,
const hp::QCollection<dim-1> &q_boundary = hp::QCollection<dim-1>(dim > 1 ?
void project (const hp::DoFHandler<dim,spacedim> &dof,
const ConstraintMatrix &constraints,
const hp::QCollection<dim> &quadrature,
- const Function<spacedim,double> &function,
+ const Function<spacedim,typename VectorType::value_type> &function,
VectorType &vec,
const bool enforce_zero_boundary = false,
const hp::QCollection<dim-1> &q_boundary = hp::QCollection<dim-1>(dim > 1 ?
template <int dim, typename VectorType, int spacedim>
void point_difference (const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
- const Function<spacedim,double> &exact_solution,
- Vector<double> &difference,
+ const Function<spacedim,typename VectorType::value_type> &exact_solution,
+ Vector<typename VectorType::value_type> &difference,
const Point<spacedim> &point);
/**
void point_difference (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
- const Function<spacedim,double> &exact_solution,
- Vector<double> &difference,
+ const Function<spacedim,typename VectorType::value_type> &exact_solution,
+ Vector<typename VectorType::value_type> &difference,
const Point<spacedim> &point);
/**
point_value (const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point,
- Vector<double> &value);
+ Vector<typename VectorType::value_type> &value);
/**
* Same as above for hp.
point_value (const hp::DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point,
- Vector<double> &value);
+ Vector<typename VectorType::value_type> &value);
/**
* Evaluate a scalar finite element function defined by the given DoFHandler
* exception of type VectorTools::ExcPointNotAvailableHere is thrown.
*/
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
point_value (const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point);
* exception of type VectorTools::ExcPointNotAvailableHere is thrown.
*/
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
point_value (const hp::DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point);
const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point,
- Vector<double> &value);
+ Vector<typename VectorType::value_type> &value);
/**
* Same as above for hp.
const hp::DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point,
- Vector<double> &value);
+ Vector<typename VectorType::value_type> &value);
/**
* Evaluate a scalar finite element function defined by the given DoFHandler
* exception of type VectorTools::ExcPointNotAvailableHere is thrown.
*/
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
point_value (const Mapping<dim,spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
* exception of type VectorTools::ExcPointNotAvailableHere is thrown.
*/
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
point_value (const hp::MappingCollection<dim,spacedim> &mapping,
const hp::DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
* mean value and subtract it right inside the evaluation routine.
*/
template <int dim, typename VectorType, int spacedim>
- double compute_mean_value (const Mapping<dim, spacedim> &mapping,
- const DoFHandler<dim,spacedim> &dof,
- const Quadrature<dim> &quadrature,
- const VectorType &v,
- const unsigned int component);
+ typename VectorType::value_type
+ compute_mean_value (const Mapping<dim, spacedim> &mapping,
+ const DoFHandler<dim,spacedim> &dof,
+ const Quadrature<dim> &quadrature,
+ const VectorType &v,
+ const unsigned int component);
/**
* Calls the other compute_mean_value() function, see above, with
* <tt>mapping=MappingQGeneric@<dim@>(1)</tt>.
*/
template <int dim, typename VectorType, int spacedim>
- double compute_mean_value (const DoFHandler<dim,spacedim> &dof,
- const Quadrature<dim> &quadrature,
- const VectorType &v,
- const unsigned int component);
+ typename VectorType::value_type
+ compute_mean_value (const DoFHandler<dim,spacedim> &dof,
+ const Quadrature<dim> &quadrature,
+ const VectorType &v,
+ const unsigned int component);
//@}
/**
* Geometrical interpolation
template <int, int> class DoFHandlerType>
void interpolate (const Mapping<dim,spacedim> &mapping,
const DoFHandlerType<dim,spacedim> &dof,
- const Function<spacedim> &function,
+ const Function<spacedim, typename VectorType::value_type> &function,
VectorType &vec)
{
+ typedef typename VectorType::value_type number;
Assert (vec.size() == dof.n_dofs(),
ExcDimensionMismatch (vec.size(), dof.n_dofs()));
Assert (dof.get_fe().n_components() == function.n_components,
// have two versions, one for system fe
// and one for scalar ones, to take the
// more efficient one respectively
- std::vector<std::vector<double> > function_values_scalar(fe.size());
- std::vector<std::vector<Vector<double> > > function_values_system(fe.size());
+ std::vector<std::vector<number> > function_values_scalar(fe.size());
+ std::vector<std::vector<Vector<number> > > function_values_system(fe.size());
// Make a quadrature rule from support points
// to feed it into FEValues
// all components
function_values_system[fe_index]
.resize (n_rep_points[fe_index],
- Vector<double> (fe[fe_index].n_components()));
+ Vector<number> (fe[fe_index].n_components()));
function.vector_value_list (rep_points,
function_values_system[fe_index]);
// distribute the function
template <typename VectorType, typename DoFHandlerType>
void interpolate (const DoFHandlerType &dof,
- const Function<DoFHandlerType::space_dimension> &function,
+ const Function<DoFHandlerType::space_dimension,typename VectorType::value_type> &function,
VectorType &vec)
{
interpolate(StaticMappingQ1<DoFHandlerType::dimension,
const InVector &data_1,
OutVector &data_2)
{
- Vector<double> cell_data_1(dof_1.get_fe().dofs_per_cell);
- Vector<double> cell_data_2(dof_2.get_fe().dofs_per_cell);
+ typedef typename OutVector::value_type number;
+ Vector<number> cell_data_1(dof_1.get_fe().dofs_per_cell);
+ Vector<number> cell_data_2(dof_2.get_fe().dofs_per_cell);
std::vector<short unsigned int> touch_count (dof_2.n_dofs(), 0); //TODO: check on datatype... kinda strange (UK)
std::vector<types::global_dof_index> local_dof_indices (dof_2.get_fe().dofs_per_cell);
std::vector<Tensor<1,spacedim> >(n_components));
}
+ namespace
+ {
+ template<typename number>
+ double mean_to_double(const number &mean_value)
+ {
+ return mean_value;
+ }
+
+ template<typename number>
+ double mean_to_double(const std::complex<number> &mean_value)
+ {
+ // we need to return double as a norm, but mean value is a complex
+ // number. Panick and return real-part while warning the user that
+ // he shall never do that.
+ Assert ( false, ExcMessage("Mean value norm is not implemented for complex-valued vectors") );
+ return mean_value.real();
+ }
+ }
+
// avoid compiling inner function for many vector types when we always
// really do the same thing by putting the main work into this helper
}
double diff = 0;
+ Number diff_mean = 0;
// First work on function values:
switch (norm)
// Compute values in quadrature points and integrate
for (unsigned int q=0; q<n_q_points; ++q)
{
- double sum = 0;
+ Number sum = 0;
for (unsigned int k=0; k<n_components; ++k)
sum += data.psi_values[q](k) * data.weight_vectors[q](k);
- diff += sum * fe_values.JxW(q);
+ diff_mean += sum * fe_values.JxW(q);
}
break;
double sum = 0;
for (unsigned int k=0; k<n_components; ++k)
sum += std::pow(
- static_cast<double>(data.psi_values[q](k)*data.psi_values[q](k)),
+ static_cast<double>(numbers::NumberTraits<Number>::abs_square(data.psi_values[q](k))),
exponent/2.) * data.weight_vectors[q](k);
diff += sum * fe_values.JxW(q);
}
{
double sum = 0;
for (unsigned int k=0; k<n_components; ++k)
- sum += data.psi_values[q](k) * data.psi_values[q](k) *
+ sum += numbers::NumberTraits<Number>::abs_square(data.psi_values[q](k)) *
data.weight_vectors[q](k);
diff += sum * fe_values.JxW(q);
}
double sum = 0;
for (unsigned int k=0; k<n_components; ++k)
sum += std::pow(
- static_cast<double>(data.psi_grads[q][k]*data.psi_grads[q][k]),
+ data.psi_grads[q][k].norm_square(),
exponent/2.) * data.weight_vectors[q](k);
diff += sum * fe_values.JxW(q);
}
{
double sum = 0;
for (unsigned int k=0; k<n_components; ++k)
- sum += (data.psi_grads[q][k] * data.psi_grads[q][k]) *
+ sum += data.psi_grads[q][k].norm_square() *
data.weight_vectors[q](k);
diff += sum * fe_values.JxW(q);
}
ExcMessage ("You can only ask for the Hdiv norm for a finite element "
"with at least 'dim' components. In that case, this function "
"will take the divergence of the first 'dim' components."));
- double sum = 0;
+ Number sum = 0;
// take the trace of the derivatives scaled by the weight and square it
for (unsigned int k=0; k<dim; ++k)
sum += data.psi_grads[q][k][k] * std::sqrt(data.weight_vectors[q](k));
- diff += sum * sum * fe_values.JxW(q);
+ diff += numbers::NumberTraits<Number>::abs_square(sum) * fe_values.JxW(q);
}
diff = std::sqrt(diff);
break;
break;
}
+ if (norm == mean)
+ diff = mean_to_double(diff_mean);
+
// append result of this cell to the end of the vector
AssertIsFinite(diff);
return diff;
void
point_difference (const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
- const Function<spacedim> &exact_function,
- Vector<double> &difference,
+ const Function<spacedim, typename VectorType::value_type> &exact_function,
+ Vector<typename VectorType::value_type> &difference,
const Point<spacedim> &point)
{
point_difference(StaticMappingQ1<dim>::mapping,
point_difference (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
- const Function<spacedim> &exact_function,
- Vector<double> &difference,
+ const Function<spacedim, typename VectorType::value_type> &exact_function,
+ Vector<typename VectorType::value_type> &difference,
const Point<spacedim> &point)
{
typedef typename VectorType::value_type Number;
point_value (const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point,
- Vector<double> &value)
+ Vector<typename VectorType::value_type> &value)
{
point_value (StaticMappingQ1<dim,spacedim>::mapping,
point_value (const hp::DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point,
- Vector<double> &value)
+ Vector<typename VectorType::value_type> &value)
{
point_value(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
dof,
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
point_value (const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point)
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
point_value (const hp::DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point)
const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point,
- Vector<double> &value)
+ Vector<typename VectorType::value_type> &value)
{
typedef typename VectorType::value_type Number;
const FiniteElement<dim> &fe = dof.get_fe();
const hp::DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
const Point<spacedim> &point,
- Vector<double> &value)
+ Vector<typename VectorType::value_type> &value)
{
typedef typename VectorType::value_type Number;
const hp::FECollection<dim, spacedim> &fe = dof.get_fe();
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
point_value (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
Assert(dof.get_fe().n_components() == 1,
ExcMessage ("Finite element is not scalar as is necessary for this function"));
- Vector<double> value(1);
+ Vector<typename VectorType::value_type> value(1);
point_value(mapping, dof, fe_function, point, value);
return value(0);
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
point_value (const hp::MappingCollection<dim, spacedim> &mapping,
const hp::DoFHandler<dim,spacedim> &dof,
const VectorType &fe_function,
Assert(dof.get_fe().n_components() == 1,
ExcMessage ("Finite element is not scalar as is necessary for this function"));
- Vector<double> value(1);
+ Vector<typename VectorType::value_type> value(1);
point_value(mapping, dof, fe_function, point, value);
return value(0);
for (unsigned int i=0; i<n; ++i)
if (p_select[i])
{
- s += v(i);
+ typename VectorType::value_type vi = v(i);
+ s += vi;
++counter;
}
// Error out if we have not constrained anything. Note that in this
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
compute_mean_value (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &quadrature,
template <int dim, typename VectorType, int spacedim>
- double
+ typename VectorType::value_type
compute_mean_value (const DoFHandler<dim,spacedim> &dof,
const Quadrature<dim> &quadrature,
const VectorType &v,
void interpolate
(const Mapping<deal_II_dimension,deal_II_space_dimension> &,
const DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
- const Function<deal_II_space_dimension>&,
+ const Function<deal_II_space_dimension, VEC::value_type>&,
VEC&);
template
void interpolate
(const DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
- const Function<deal_II_space_dimension>&,
+ const Function<deal_II_space_dimension, VEC::value_type>&,
VEC&);
template
void interpolate
(const Mapping<deal_II_dimension>&,
const hp::DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
+ const Function<deal_II_dimension, VEC::value_type>&,
VEC&);
template
void interpolate
(const hp::DoFHandler<deal_II_dimension>&,
- const Function<deal_II_dimension>&,
+ const Function<deal_II_dimension, VEC::value_type>&,
VEC&);
template
namespace VectorTools \{
template
- double compute_mean_value<deal_II_dimension>
+ VEC::value_type compute_mean_value<deal_II_dimension>
(const Mapping<deal_II_dimension,deal_II_space_dimension>&,
const DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
const Quadrature<deal_II_dimension>&,
const unsigned int);
template
- double compute_mean_value<deal_II_dimension>
+ VEC::value_type compute_mean_value<deal_II_dimension>
(const DoFHandler<deal_II_dimension,deal_II_space_dimension>&,
const Quadrature<deal_II_dimension>&,
const VEC&,
const hp::DoFHandler<deal_II_dimension>&,
const VEC&,
const Point<deal_II_dimension>&,
- Vector<double>&);
+ Vector<VEC::value_type>&);
template
- double point_value<deal_II_dimension> (
+ VEC::value_type point_value<deal_II_dimension> (
const hp::DoFHandler<deal_II_dimension>&,
const VEC&,
const Point<deal_II_dimension>&);
const hp::DoFHandler<deal_II_dimension>&,
const VEC&,
const Point<deal_II_dimension>&,
- Vector<double>&);
+ Vector<VEC::value_type>&);
template
- double point_value<deal_II_dimension> (
+ VEC::value_type point_value<deal_II_dimension> (
const hp::MappingCollection<deal_II_dimension>&,
const hp::DoFHandler<deal_II_dimension>&,
const VEC&,
void point_difference<deal_II_dimension> (
const DoFHandler<deal_II_dimension>&,
const VEC&,
- const Function<deal_II_dimension>&,
- Vector<double>&,
+ const Function<deal_II_dimension, VEC::value_type>&,
+ Vector<VEC::value_type>&,
const Point<deal_II_dimension>&);
template
const Mapping<deal_II_dimension>&,
const DoFHandler<deal_II_dimension>&,
const VEC&,
- const Function<deal_II_dimension>&,
- Vector<double>&,
+ const Function<deal_II_dimension,VEC::value_type>&,
+ Vector<VEC::value_type>&,
const Point<deal_II_dimension>&);
template
const DoFHandler<deal_II_dimension>&,
const VEC&,
const Point<deal_II_dimension>&,
- Vector<double>&);
+ Vector<VEC::value_type>&);
template
- double point_value<deal_II_dimension> (
+ VEC::value_type point_value<deal_II_dimension> (
const DoFHandler<deal_II_dimension>&,
const VEC&,
const Point<deal_II_dimension>&);
const DoFHandler<deal_II_dimension>&,
const VEC&,
const Point<deal_II_dimension>&,
- Vector<double>&);
+ Vector<VEC::value_type>&);
template
- double point_value<deal_II_dimension> (
+ VEC::value_type point_value<deal_II_dimension> (
const Mapping<deal_II_dimension>&,
const DoFHandler<deal_II_dimension>&,
const VEC&,