point_gradient (const DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point,
- std::vector<Tensor<1, spacedim> > &value);
+ std::vector<Tensor<1, spacedim, typename InVector::value_type> > &value);
/**
* Same as above for hp.
point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point,
- std::vector<Tensor<1, spacedim> > &value);
+ std::vector<Tensor<1, spacedim, typename InVector::value_type> > &value);
/**
* Evaluate a scalar finite element function defined by the given DoFHandler
* exception of type VectorTools::ExcPointNotAvailableHere is thrown.
*/
template <int dim, class InVector, int spacedim>
- Tensor<1, spacedim>
+ Tensor<1, spacedim, typename InVector::value_type>
point_gradient (const DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point);
* exception of type VectorTools::ExcPointNotAvailableHere is thrown.
*/
template <int dim, class InVector, int spacedim>
- Tensor<1, spacedim>
+ Tensor<1, spacedim, typename InVector::value_type>
point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point);
const DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point,
- std::vector<Tensor<1, spacedim> > &value);
+ std::vector<Tensor<1, spacedim, typename InVector::value_type> > &value);
/**
* Same as above for hp.
const hp::DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point,
- std::vector<Tensor<1, spacedim> > &value);
+ std::vector<Tensor<1, spacedim, typename InVector::value_type> > &value);
/**
* Evaluate a scalar finite element function defined by the given DoFHandler
* exception of type VectorTools::ExcPointNotAvailableHere is thrown.
*/
template <int dim, class InVector, int spacedim>
- Tensor<1, spacedim>
+ Tensor<1, spacedim, typename InVector::value_type>
point_gradient (const Mapping<dim,spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
* exception of type VectorTools::ExcPointNotAvailableHere is thrown.
*/
template <int dim, class InVector, int spacedim>
- Tensor<1, spacedim>
+ Tensor<1, spacedim, typename InVector::value_type>
point_gradient (const hp::MappingCollection<dim,spacedim> &mapping,
const hp::DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
namespace internal
{
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
struct IDScratchData
{
IDScratchData (const dealii::hp::MappingCollection<dim,spacedim> &mapping,
void resize_vectors (const unsigned int n_q_points,
const unsigned int n_components);
- std::vector<dealii::Vector<double> > function_values;
- std::vector<std::vector<Tensor<1,spacedim> > > function_grads;
+ std::vector<dealii::Vector<Number> > function_values;
+ std::vector<std::vector<Tensor<1,spacedim,Number> > > function_grads;
std::vector<double> weight_values;
std::vector<dealii::Vector<double> > weight_vectors;
- std::vector<dealii::Vector<double> > psi_values;
- std::vector<std::vector<Tensor<1,spacedim> > > psi_grads;
- std::vector<double> psi_scalar;
+ std::vector<dealii::Vector<Number> > psi_values;
+ std::vector<std::vector<Tensor<1,spacedim,Number> > > psi_grads;
+ std::vector<Number> psi_scalar;
std::vector<double> tmp_values;
+ std::vector<dealii::Vector<double> > tmp_vector_values;
std::vector<Tensor<1,spacedim> > tmp_gradients;
+ std::vector<std::vector<Tensor<1,spacedim> > > tmp_vector_gradients;
dealii::hp::FEValues<dim,spacedim> x_fe_values;
};
- template <int dim, int spacedim>
- IDScratchData<dim,spacedim>
+ template <int dim, int spacedim, typename Number>
+ IDScratchData<dim,spacedim,Number>
::IDScratchData(const dealii::hp::MappingCollection<dim,spacedim> &mapping,
const dealii::hp::FECollection<dim,spacedim> &fe,
const dealii::hp::QCollection<dim> &q,
x_fe_values(mapping, fe, q, update_flags)
{}
- template <int dim, int spacedim>
- IDScratchData<dim,spacedim>::IDScratchData (const IDScratchData &data)
+ template <int dim, int spacedim, typename Number>
+ IDScratchData<dim,spacedim,Number>::IDScratchData (const IDScratchData &data)
:
x_fe_values(data.x_fe_values.get_mapping_collection(),
data.x_fe_values.get_fe_collection(),
data.x_fe_values.get_update_flags())
{}
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
void
- IDScratchData<dim,spacedim>::resize_vectors (const unsigned int n_q_points,
- const unsigned int n_components)
+ IDScratchData<dim,spacedim,Number>::resize_vectors (const unsigned int n_q_points,
+ const unsigned int n_components)
{
function_values.resize (n_q_points,
- dealii::Vector<double>(n_components));
+ dealii::Vector<Number>(n_components));
function_grads.resize (n_q_points,
- std::vector<Tensor<1,spacedim> >(n_components));
+ std::vector<Tensor<1,spacedim,Number> >(n_components));
weight_values.resize (n_q_points);
weight_vectors.resize (n_q_points,
dealii::Vector<double>(n_components));
psi_values.resize (n_q_points,
- dealii::Vector<double>(n_components));
+ dealii::Vector<Number>(n_components));
psi_grads.resize (n_q_points,
- std::vector<Tensor<1,spacedim> >(n_components));
+ std::vector<Tensor<1,spacedim,Number> >(n_components));
psi_scalar.resize (n_q_points);
tmp_values.resize (n_q_points);
+ tmp_vector_values.resize (n_q_points,
+ dealii::Vector<double>(n_components));
tmp_gradients.resize (n_q_points);
+ tmp_vector_gradients.resize (n_q_points,
+ std::vector<Tensor<1,spacedim> >(n_components));
}
// avoid compiling inner function for many vector types when we always
// really do the same thing by putting the main work into this helper
// function
- template <int dim, int spacedim>
+ template <int dim, int spacedim, typename Number>
double
integrate_difference_inner (const Function<spacedim> &exact_solution,
const NormType &norm,
const UpdateFlags update_flags,
const double exponent,
const unsigned int n_components,
- IDScratchData<dim,spacedim> &data)
+ IDScratchData<dim,spacedim,Number> &data)
{
const bool fe_is_system = (n_components != 1);
const dealii::FEValues<dim, spacedim> &fe_values = data.x_fe_values.get_present_fe_values ();
if (update_flags & update_values)
{
// first compute the exact solution (vectors) at the quadrature
- // points try to do this as efficient as possible by avoiding a
+ // points. try to do this as efficient as possible by avoiding a
// second virtual function call in case the function really has only
// one component
+ //
+ // TODO: we have to work a bit here because the Function<dim,double>
+ // interface of the argument denoting the exact function only
+ // provides us with double/Tensor<1,dim> values, rather than
+ // with the correct data type. so evaluate into a temp
+ // object, then copy around
if (fe_is_system)
- exact_solution.vector_value_list (fe_values.get_quadrature_points(),
- data.psi_values);
+ {
+ exact_solution.vector_value_list (fe_values.get_quadrature_points(),
+ data.tmp_vector_values);
+ for (unsigned int i=0; i<n_q_points; ++i)
+ data.psi_values[i] = data.tmp_vector_values[i];
+ }
else
{
exact_solution.value_list (fe_values.get_quadrature_points(),
// then subtract finite element fe_function
for (unsigned int q=0; q<n_q_points; ++q)
- data.psi_values[q] -= data.function_values[q];
+ for (unsigned int i=0; i<data.psi_values[q].size(); ++i)
+ data.psi_values[q][i] -= data.function_values[q][i];
}
// Do the same for gradients, if required
// calls when calling gradient_list for functions that are really
// scalar functions
if (fe_is_system)
- exact_solution.vector_gradient_list (fe_values.get_quadrature_points(),
- data.psi_grads);
+ {
+ exact_solution.vector_gradient_list (fe_values.get_quadrature_points(),
+ data.tmp_vector_gradients);
+ for (unsigned int i=0; i<n_q_points; ++i)
+ for (unsigned int comp=0; comp<data.psi_grads[i].size(); ++comp)
+ data.psi_grads[i][comp] = data.tmp_vector_gradients[i][comp];
+ }
else
{
exact_solution.gradient_list (fe_values.get_quadrature_points(),
if (update_flags & update_normal_vectors)
for (unsigned int k=0; k<n_components; ++k)
for (unsigned int q=0; q<n_q_points; ++q)
- data.psi_grads[q][k] -= (data.function_grads[q][k] +
- (data.psi_grads[q][k] * // (f.n) n
- Tensor<1,spacedim>(fe_values.normal_vector(q)))*
- fe_values.normal_vector(q));
+ {
+ // compute (f.n) n
+ const Number f_dot_n
+ = (data.psi_grads[q][k] * Tensor<1,spacedim,Number>(fe_values.normal_vector(q)));
+ const Tensor<1,spacedim,Number> f_dot_n_times_n
+ = f_dot_n * Tensor<1,spacedim,Number>(fe_values.normal_vector(q));
+
+ data.psi_grads[q][k] -= (data.function_grads[q][k] + f_dot_n_times_n);
+ }
else
for (unsigned int k=0; k<n_components; ++k)
for (unsigned int q=0; q<n_q_points; ++q)
- data.psi_grads[q][k] -= data.function_grads[q][k];
+ for (unsigned int d=0; d<spacedim; ++d)
+ data.psi_grads[q][k][d] -= data.function_grads[q][k][d];
}
double diff = 0;
case W1infty_norm:
for (unsigned int q=0; q<n_q_points; ++q)
for (unsigned int k=0; k<n_components; ++k)
- diff = std::max (diff, std::abs(data.psi_values[q](k)*
- data.weight_vectors[q](k)));
+ diff = std::max (diff, double(std::abs(data.psi_values[q](k)*
+ data.weight_vectors[q](k))));
break;
case H1_seminorm:
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));
+ t = std::max(t,
+ double(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;
}
dealii::hp::FECollection<dim,spacedim> fe_collection (dof.get_fe());
- IDScratchData<dim,spacedim> data(mapping, fe_collection, q, update_flags);
-
- // FIXME
- // temporary vectors of consistent with InVector type.
- // Need these because IDScratchData does not have a template type for the InVector
- std::vector<dealii::Vector<Number>> function_values;
- std::vector<std::vector<Tensor<1,spacedim,Number> > > function_grads;
+ IDScratchData<dim,spacedim,typename InVector::value_type> data(mapping, fe_collection, q, update_flags);
// loop over all cells
for (typename DH::active_cell_iterator cell = dof.begin_active();
const unsigned int n_q_points = fe_values.n_quadrature_points;
data.resize_vectors (n_q_points, n_components);
- //FIXME:
- function_values.resize (n_q_points,
- dealii::Vector<Number>(n_components));
- function_grads.resize (n_q_points,
- std::vector<Tensor<1,spacedim,Number> >(n_components));
-
if (update_flags & update_values)
- fe_values.get_function_values (fe_function, function_values);
+ fe_values.get_function_values (fe_function, data.function_values);
if (update_flags & update_gradients)
- fe_values.get_function_gradients (fe_function, function_grads);
-
- // FIXME
- for (unsigned int q = 0; q < n_q_points; q++)
- for (unsigned int c = 0; c < n_components; c++)
- {
- data.function_values[q][c] = function_values[q][c];
- data.function_grads[q][c] = function_grads[q][c];
- }
+ fe_values.get_function_gradients (fe_function, data.function_grads);
difference(cell->active_cell_index()) =
- integrate_difference_inner (exact_solution, norm, weight,
- update_flags, exponent,
- n_components, data);
+ integrate_difference_inner<dim,spacedim,typename InVector::value_type> (exact_solution, norm, weight,
+ update_flags, exponent,
+ n_components, data);
}
else
// the cell is a ghost cell or is artificial. write a zero into the
point_gradient (const DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point,
- std::vector<Tensor<1, spacedim> > &gradients)
+ std::vector<Tensor<1, spacedim, typename InVector::value_type> > &gradients)
{
point_gradient (StaticMappingQ1<dim,spacedim>::mapping,
point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point,
- std::vector<Tensor<1, spacedim> > &gradients)
+ std::vector<Tensor<1, spacedim, typename InVector::value_type> > &gradients)
{
point_gradient(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
dof,
template <int dim, class InVector, int spacedim>
- Tensor<1, spacedim>
+ Tensor<1, spacedim, typename InVector::value_type>
point_gradient (const DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point)
template <int dim, class InVector, int spacedim>
- Tensor<1, spacedim>
+ Tensor<1, spacedim, typename InVector::value_type>
point_gradient (const hp::DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point)
const DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point,
- std::vector<Tensor<1, spacedim> > &gradient)
+ std::vector<Tensor<1, spacedim, typename InVector::value_type> > &gradient)
{
const FiniteElement<dim> &fe = dof.get_fe();
// the given fe_function at this point
typedef typename InVector::value_type Number;
std::vector<std::vector<Tensor<1, dim, Number> > >
- u_gradient(1, std::vector<Tensor<1, dim, Number> > (fe.n_components()));
+ u_gradient(1, std::vector<Tensor<1, dim, Number> > (fe.n_components()));
fe_values.get_function_gradients(fe_function, u_gradient);
gradient = u_gradient[0];
const hp::DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
const Point<spacedim> &point,
- std::vector<Tensor<1, spacedim> > &gradient)
+ std::vector<Tensor<1, spacedim, typename InVector::value_type> > &gradient)
{
typedef typename InVector::value_type Number;
const hp::FECollection<dim, spacedim> &fe = dof.get_fe();
// the given fe_function at this point
typedef typename InVector::value_type Number;
std::vector<std::vector<Tensor<1, dim, Number> > >
- u_gradient(1, std::vector<Tensor<1, dim, Number> > (fe.n_components()));
+ u_gradient(1, std::vector<Tensor<1, dim, Number> > (fe.n_components()));
fe_values.get_function_gradients(fe_function, u_gradient);
gradient = u_gradient[0];
template <int dim, class InVector, int spacedim>
- Tensor<1, spacedim>
+ Tensor<1, spacedim, typename InVector::value_type>
point_gradient (const Mapping<dim, spacedim> &mapping,
const DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
Assert(dof.get_fe().n_components() == 1,
ExcMessage ("Finite element is not scalar as is necessary for this function"));
- std::vector<Tensor<1, dim> > gradient(1);
+ std::vector<Tensor<1, dim, typename InVector::value_type> > gradient(1);
point_gradient (mapping, dof, fe_function, point, gradient);
return gradient[0];
}
+
template <int dim, class InVector, int spacedim>
- Tensor<1, spacedim>
+ Tensor<1, spacedim, typename InVector::value_type>
point_gradient (const hp::MappingCollection<dim, spacedim> &mapping,
const hp::DoFHandler<dim,spacedim> &dof,
const InVector &fe_function,
Assert(dof.get_fe().n_components() == 1,
ExcMessage ("Finite element is not scalar as is necessary for this function"));
- std::vector<Tensor<1, dim> > gradient(1);
+ std::vector<Tensor<1, dim, typename InVector::value_type> > gradient(1);
point_gradient (mapping, dof, fe_function, point, gradient);
return gradient[0];
namespace internal
{
template <typename Number>
- void set_possibly_complex_number(Number &n, const double &r, const double &i)
+ void set_possibly_complex_number(Number &n, const double &r, const double &)
{
n = r;
}