* usual in C++, the rightmost
* index marches fastest.
*/
- void unroll(vector<double> & result) const;
+ void unroll(Vector<double> & result) const;
/**
*/
void clear ();
- /**
- * Exception
- */
- DeclException2(ExcWrongVectorSize, unsigned, int, << "Tensor has " << arg1
- << " entries, but vector has size " << arg2);
+// /**
+// * Exception
+// */
+// DeclException2(ExcWrongVectorSize, unsigned, int, << "Tensor has " << arg1
+// << " entries, but vector has size " << arg2);
- /**
- * Exception
- */
- DeclException1 (ExcInvalidIndex,
- int,
- << "Invalid index " << arg1);
+// /**
+// * Exception
+// */
+// DeclException1 (ExcInvalidIndex,
+// int,
+// << "Invalid index " << arg1);
private:
/**
* Array of tensors holding the
/**
* Help function for unroll.
*/
- void unroll_recursion(vector<double> & result, unsigned& start_index) const;
+ void unroll_recursion(Vector<double> & result, unsigned& start_index) const;
template<>
- friend void Tensor<rank_+1,dim>::unroll_recursion(vector<double> &, unsigned&) const;
+ friend void Tensor<rank_+1,dim>::unroll_recursion(Vector<double> &, unsigned&) const;
};
#include <base/exceptions.h>
#include <iostream>
#include <vector>
-
-
+#include <lac/forward-declarations.h>
// general template; specialized for rank==1; the general template is in
* usual in C++, the rightmost
* index marches fastest.
*/
- void unroll(vector<double> & result) const;
+ void unroll(Vector<double> & result) const;
protected:
/**
/**
* Help function for unroll.
*/
- void unroll_recursion(vector<double> & result, unsigned& start_index) const;
+ void unroll_recursion(Vector<double> & result, unsigned& start_index) const;
template<>
- friend void Tensor<2,dim>::unroll_recursion(vector<double> &, unsigned&) const;
+ friend void Tensor<2,dim>::unroll_recursion(Vector<double> &, unsigned&) const;
};
-
-
-
+ /**
+ * Exception
+ */
+DeclException2(ExcWrongVectorSize, unsigned, int, << "Tensor has " << arg1
+ << " entries, but vector has size " << arg2);
DeclException1 (ExcDimTooSmall,
int,
<< "Given dimensions must be >= 1, but was " << arg1);
#include <base/point.h>
#include <base/functiontime.h>
#include <base/tensorindex.h>
-
+#include <lac/forward-declarations.h>
template<int rank_, int dim>
class Tensor;
* the same size as the #n_components#
* array.
*/
- virtual void value (const Point<dim> &p, vector<double> &values) const;
+ virtual void value (const Point<dim> &p, Vector<double> &values) const;
/**
* Set #values# to the point values
* array.
*/
virtual void value_list (const vector<Point<dim> > &points,
- vector<vector<double> > &values) const;
+ vector<Vector<double> > &values) const;
/**
* Set #gradients# to the gradients of
virtual void gradient_list (const vector<Point<dim> > &points,
vector<Tensor<rank_+1,dim> > &gradients) const;
+ /**
+ * See #VectorFunction#.
+ */
+ virtual void value (const Point<dim> &points,
+ Vector<double> &values) const;
+
/**
* See #VectorFunction#.
*/
virtual void value_list (const vector<Point<dim> > &points,
- vector<vector<double> > &values) const;
+ vector<Vector<double> > &values) const;
/**
* See #VectorFunction#.
#include <base/tensor.h>
#include <cmath>
+#include <lac/vector.h>
-template <int dim> void
-Tensor<1,dim>::unroll( vector<double>& result) const
+template <int dim>
+void
+Tensor<1,dim>::unroll( Vector<double>& result) const
{
- Assert(false, ExcNotImplemented());
+ Assert(result.size()==dim,
+ ExcWrongVectorSize(dim, result.size()));
+
+ unsigned index = 0;
+ unroll_recursion(result,index);
}
-template <int rank_, int dim> void
-Tensor<rank_, dim>::unroll( vector<double>& result) const
+template <int rank_, int dim>
+void
+Tensor<rank_, dim>::unroll( Vector<double>& result) const
{
Assert(result.size()==pow(dim,rank_),
ExcWrongVectorSize(pow(dim,rank_), result.size()));
}
-template <int rank_, int dim> void
-Tensor<rank_, dim>::unroll_recursion( vector<double>& result, unsigned& index) const
+template <int rank_, int dim>
+void
+Tensor<rank_, dim>::unroll_recursion( Vector<double>& result, unsigned& index) const
{
for (unsigned i=0; i<dim; ++i)
{
}
-template<int dim> void
-Tensor<1,dim>::unroll_recursion( vector<double>& result, unsigned& index) const
+template<int dim>
+void
+Tensor<1,dim>::unroll_recursion( Vector<double>& result, unsigned& index) const
{
for (unsigned i=0; i<dim; ++i)
{
- cerr << "[" << index << ',' << operator[](i) << ']' << endl;
- result[index++] = operator[](i);
+// cerr << "[" << index << ',' << operator[](i) << ']' << endl;
+ result(index++) = operator[](i);
}
}
#include <vector>
#include <base/tensor.h>
#include <cmath>
+#include <lac/vector.h>
template <int dim>
VectorFunction<dim>::VectorFunction(unsigned n_components, const double initial_time)
*/
template <int dim> void
-VectorFunction<dim>::value (const Point<dim> &, vector<double> &) const
+VectorFunction<dim>::value (const Point<dim> &, Vector<double> &) const
{
Assert (false, ExcPureFunctionCalled());
}
template <int dim> void
VectorFunction<dim>::value_list (const vector<Point<dim> > &,
- vector<vector<double> > &) const
+ vector<Vector<double> > &) const
{
Assert (false, ExcPureFunctionCalled());
}
template <int rank_, int dim> void
TensorFunction<rank_, dim>::value_list (const vector<Point<dim> > & points,
- vector<vector<double> > & values) const
+ vector<Vector<double> > & values) const
{
Assert (values.size() == points.size(),
ExcVectorHasWrongSize(values.size(), points.size()));
* elements.
*/
void get_function_values (const Vector<double> &fe_function,
- vector<vector<double> > &values) const;
+ vector<Vector<double> > &values) const;
/**
* Return the gradient of the #i#th shape
* much needed here and saving memory may
* be a good goal when using many cells.
*/
- void refine (const Vector<float> &criteria,
+ template <typename number>
+ void refine (const Vector<number> &criteria,
const double threshold);
/**
* much needed here and saving memory may
* be a good goal when using many cells.
*/
- void coarsen (const Vector<float> &criteria,
+ template <typename number>
+ void coarsen (const Vector<number> &criteria,
const double threshold);
/**
* much needed here and saving memory may
* be a good goal when using many cells.
*/
- void refine_and_coarsen_fixed_number (const Vector<float> &criteria,
+ template <typename number>
+ void refine_and_coarsen_fixed_number (const Vector<number> &criteria,
const double top_fraction_of_cells,
const double bottom_fraction_of_cells);
* much needed here and saving memory may
* be a good goal when using many cells.
*/
- void refine_and_coarsen_fixed_fraction (const Vector<float> &criteria,
+ template<typename number>
+ void refine_and_coarsen_fixed_fraction (const Vector<number> &criteria,
const double top_fraction,
const double bottom_fraction);
* Compute the error for the solution of a system.
* See the other #integrate_difference#.
*/
-// static void integrate_difference (const DoFHandler<dim> &dof,
-// const Vector<double> &fe_function,
-// const TensorFunction<1,dim>&exact_solution,
-// Vector<float> &difference,
-// const Quadrature<dim> &q,
-// const FiniteElement<dim> &fe,
-// const NormType &norm,
-// const Boundary<dim> &boundary);
+ static void integrate_difference (const DoFHandler<dim> &dof,
+ const Vector<double> &fe_function,
+ const VectorFunction<dim>&exact_solution,
+ Vector<float> &difference,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe,
+ const NormType &norm,
+ const Boundary<dim> &boundary);
/**
// 1. Vertices
total_index = 0;
- for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<2>::vertices_per_cell ;
+ for (unsigned vertex_number= 0 ; vertex_number < GeometryInfo<dim>::vertices_per_cell ;
++vertex_number)
{
unsigned comp_start = 0;
}
// 2. Lines
- for (unsigned line_number= 0 ; line_number < GeometryInfo<2>::lines_per_cell ;
+ for (unsigned line_number= 0 ; line_number < GeometryInfo<dim>::lines_per_cell ;
++line_number)
{
unsigned comp_start = 0;
}
// 3. Quads
- for (unsigned quad_number= 0 ; quad_number < GeometryInfo<2>::quads_per_cell ;
+ for (unsigned quad_number= 0 ; quad_number < GeometryInfo<dim>::quads_per_cell ;
++quad_number)
{
unsigned comp_start = 0;
}
// 4. Hex
- for (unsigned hex_number= 0 ; hex_number < GeometryInfo<2>::hexes_per_cell ;
+ for (unsigned hex_number= 0 ; hex_number < GeometryInfo<dim>::hexes_per_cell ;
++hex_number)
{
unsigned comp_start = 0;
template <int dim>
void FEValuesBase<dim>::get_function_values (const Vector<double> &fe_function,
- vector< vector<double> > &values) const
+ vector< Vector<double> > &values) const
{
Assert (fe->n_components == values.size(),
ExcWrongNoOfComponents());
// functions
for (unsigned int point=0; point<n_quadrature_points; ++point)
for (unsigned int shape_func=0; shape_func<total_dofs; ++shape_func)
- values[fe->system_to_component_index(shape_func).first][point]
+ values[point](fe->system_to_component_index(shape_func).first)
+= (dof_values(shape_func) * shape_values[selected_dataset](shape_func, point));
};
-template <int dim>
-void Triangulation<dim>::refine (const Vector<float> &criteria,
- const double threshold) {
+template <int dim> template <typename number>
+void Triangulation<dim>::refine (const Vector<number> &criteria,
+ const double threshold)
+{
Assert (criteria.size() == n_active_cells(),
ExcInvalidVectorSize(criteria.size(), n_active_cells()));
-template <int dim>
-void Triangulation<dim>::coarsen (const Vector<float> &criteria,
- const double threshold) {
+template <int dim> template <typename number>
+void Triangulation<dim>::coarsen (const Vector<number> &criteria,
+ const double threshold)
+{
Assert (criteria.size() == n_active_cells(),
ExcInvalidVectorSize(criteria.size(), n_active_cells()));
-template <int dim>
-void Triangulation<dim>::refine_and_coarsen_fixed_number (const Vector<float> &criteria,
- const double top_fraction,
- const double bottom_fraction) {
+template <int dim> template <typename number>
+void
+Triangulation<dim>::refine_and_coarsen_fixed_number (const Vector<number> &criteria,
+ const double top_fraction,
+ const double bottom_fraction)
+{
// correct number of cells is
// checked in #refine#
Assert ((top_fraction>=0) && (top_fraction<=1), ExcInvalidParameterValue());
const int coarsen_cells = max(static_cast<int>(bottom_fraction*criteria.size()),
1);
- Vector<float> tmp(criteria);
+ Vector<number> tmp(criteria);
nth_element (tmp.begin(), tmp.begin()+refine_cells,
tmp.end(),
greater<double>());
-template <int dim>
+template <int dim> template <typename number>
void
-Triangulation<dim>::refine_and_coarsen_fixed_fraction (const Vector<float> &criteria,
+Triangulation<dim>::refine_and_coarsen_fixed_fraction (const Vector<number> &criteria,
const double top_fraction,
const double bottom_fraction) {
// correct number of cells is
// error, which is what we have to sum
// up and compare with
// #fraction_of_error*total_error#.
- Vector<float> tmp(criteria);
+ Vector<number> tmp(criteria);
const double total_error = tmp.l1_norm();
Vector<float> partial_sums(criteria.size());
// explicit instantiations
template class Triangulation<deal_II_dimension>;
+template void Triangulation<deal_II_dimension>
+::refine (const Vector<float> &, const double);
+
+template void Triangulation<deal_II_dimension>
+::refine (const Vector<double> &, const double);
+
+template void Triangulation<deal_II_dimension>
+::coarsen (const Vector<float> &, const double);
+
+template void Triangulation<deal_II_dimension>
+::coarsen (const Vector<double> &, const double);
+
+
+template void Triangulation<deal_II_dimension>
+::refine_and_coarsen_fixed_number (const Vector<double> &,
+ const double top_fraction,
+ const double bottom_fraction);
+
+template void Triangulation<deal_II_dimension>
+::refine_and_coarsen_fixed_number (const Vector<float> &criteria,
+ const double top_fraction,
+ const double bottom_fraction);
+
+template void Triangulation<deal_II_dimension>
+::refine_and_coarsen_fixed_fraction (const Vector<double> &criteria,
+ const double top_fraction,
+ const double bottom_fraction);
+
+template void Triangulation<deal_II_dimension>
+::refine_and_coarsen_fixed_fraction (const Vector<float> &criteria,
+ const double top_fraction,
+ const double bottom_fraction);
+
FEValues<dim> fe(dofs->get_fe(), points, UpdateFlags(update_q_points));
const StraightBoundary<dim> boundary;
- vector< vector <vector<double> > >
- values (data.size(), vector< vector<double> >(dofs->get_fe().n_components, vector<double>(points.n_quadrature_points)));
+ vector< vector <Vector<double> > >
+ values (data.size(),
+ vector< Vector<double> >(points.n_quadrature_points,
+ Vector<double>(dofs->get_fe().n_components
+ )));
for (cell=dofs->begin_active(); cell!=endc; ++cell)
{
out << pt << " ";
for (unsigned int i=0; i!=data.size(); ++i)
for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
- out << values[i][j][supp_pt]
+ out << values[i][supp_pt](j)
<< ' ';
out << endl;
};
for (unsigned int i=0; i!=data.size(); ++i)
for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
- out << values[i][j][supp_pt]
+ out << values[i][supp_pt](j)
<< ' ';
out << endl;
}
for (unsigned int i=0; i!=data.size(); ++i)
for (unsigned int j=0; j < dofs->get_fe().n_components; ++j)
- out << values[i][j][supp_pt]
+ out << values[i][supp_pt](j)
<< ' ';
out << endl;
}
#include <base/function.h>
+#include <base/tensorfunction.h>
#include <grid/dof.h>
#include <grid/dof_accessor.h>
#include <grid/tria_iterator.h>
-// template <int dim>
-// void VectorTools<dim>::integrate_difference (const DoFHandler<dim> &dof,
-// const Vector<double> &fe_function,
-// const TensorFunction<1,dim>&exact_solution,
-// Vector<float> &difference,
-// const Quadrature<dim> &q,
-// const FiniteElement<dim> &fe,
-// const NormType &norm,
-// const Boundary<dim> &boundary) {
-// Assert (fe == dof.get_fe(), ExcInvalidFE());
+template <int dim>
+void
+VectorTools<dim>::integrate_difference (const DoFHandler<dim> &dof,
+ const Vector<double> &fe_function,
+ const VectorFunction<dim>&exact_solution,
+ Vector<float> &difference,
+ const Quadrature<dim> &q,
+ const FiniteElement<dim> &fe,
+ const NormType &norm,
+ const Boundary<dim> &boundary)
+{
+ Assert (fe == dof.get_fe(), ExcInvalidFE());
-// difference.reinit (dof.get_tria().n_active_cells());
+ difference.reinit (dof.get_tria().n_active_cells());
-// UpdateFlags update_flags = UpdateFlags (update_q_points |
-// update_JxW_values);
-// if ((norm==H1_seminorm) || (norm==H1_norm))
-// update_flags = UpdateFlags (update_flags | update_gradients);
-// FEValues<dim> fe_values(fe, q, update_flags);
+ UpdateFlags update_flags = UpdateFlags (update_q_points |
+ update_JxW_values);
+ if ((norm==H1_seminorm) || (norm==H1_norm))
+ update_flags = UpdateFlags (update_flags | update_gradients);
+ FEValues<dim> fe_values(fe, q, update_flags);
-// // loop over all cells
-// DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
-// endc = dof.end();
-// for (unsigned int index=0; cell != endc; ++cell, ++index)
-// {
-// double diff=0;
-// // initialize for this cell
-// fe_values.reinit (cell, boundary);
-
-// switch (norm)
-// {
-// case mean:
-// case L1_norm:
-// case L2_norm:
-// case Linfty_norm:
-// case H1_norm:
-// {
-// // we need the finite element
-// // function \psi at the different
-// // integration points. Compute
-// // it like this:
-// // \psi(x_j)=\sum_i v_i \phi_i(x_j)
-// // with v_i the nodal values of the
-// // fe_function and \phi_i(x_j) the
-// // matrix of the trial function
-// // values at the integration point
-// // x_j. Then the vector
-// // of the \psi(x_j) is v*Phi with
-// // v being the vector of nodal
-// // values on this cell and Phi
-// // the matrix.
-// //
-// // we then need the difference:
-// // reference_function(x_j)-\psi_j
-// // and assign that to the vector
-// // \psi.
-// const unsigned int n_q_points = q.n_quadrature_points;
-// vector<double> psi (n_q_points);
-
-// // in praxi: first compute
-// // exact fe_function vector
-// exact_solution.value_list (fe_values.get_quadrature_points(),
-// psi);
-// // then subtract finite element
-// // fe_function
-// if (true)
-// {
-// vector< vector<double> > function_values (fe.n_components, vector<double>(n_q_points, 0));
-// fe_values.get_function_values (fe_function, function_values);
-
-// transform (psi.begin(), psi.end(),
-// function_values.begin(),
-// psi.begin(),
-// minus<double>());
-// };
-
-// // for L1_norm and Linfty_norm:
-// // take absolute
-// // value, for the L2_norm take
-// // square of psi
-// switch (norm)
-// {
-// case mean:
-// break;
-// case L1_norm:
-// case Linfty_norm:
-// transform (psi.begin(), psi.end(),
-// psi.begin(), ptr_fun(fabs));
-// break;
-// case L2_norm:
-// case H1_norm:
-// transform (psi.begin(), psi.end(),
-// psi.begin(), ptr_fun(sqr));
-// break;
-// default:
-// Assert (false, ExcNotImplemented());
-// };
-
-// // ok, now we have the integrand,
-// // let's compute the integral,
-// // which is
-// // sum_j psi_j JxW_j
-// // (or |psi_j| or |psi_j|^2
-// switch (norm)
-// {
-// case mean:
-// case L1_norm:
-// diff = inner_product (psi.begin(), psi.end(),
-// fe_values.get_JxW_values().begin(),
-// 0.0);
-// break;
-// case L2_norm:
-// case H1_norm:
-// diff = sqrt(inner_product (psi.begin(), psi.end(),
-// fe_values.get_JxW_values().begin(),
-// 0.0));
-// break;
-// case Linfty_norm:
-// diff = *max_element (psi.begin(), psi.end());
-// break;
-// default:
-// Assert (false, ExcNotImplemented());
-// };
-
-// // note: the H1_norm uses the result
-// // of the L2_norm and control goes
-// // over to the next case statement!
-// if (norm != H1_norm)
-// break;
-// };
-
-// case H1_seminorm:
-// {
-// // note: the computation of the
-// // H1_norm starts at the previous
-// // case statement, but continues
-// // here!
-
-// // for H1_norm: re-square L2_norm.
-// diff = sqr(diff);
-
-// // same procedure as above, but now
-// // psi is a vector of gradients
-// const unsigned int n_q_points = q.n_quadrature_points;
-// vector<Tensor<1,dim> > psi (n_q_points);
-
-// // in praxi: first compute
-// // exact fe_function vector
-// exact_solution.gradient_list (fe_values.get_quadrature_points(),
-// psi);
-
-// // then subtract finite element
-// // fe_function
-// if (true)
-// {
-// vector<Tensor<1,dim> > function_grads (n_q_points, Tensor<1,dim>());
-// fe_values.get_function_grads (fe_function, function_grads);
-
-// transform (psi.begin(), psi.end(),
-// function_grads.begin(),
-// psi.begin(),
-// minus<Tensor<1,dim> >());
-// };
-// // take square of integrand
-// vector<double> psi_square (psi.size(), 0.0);
-// for (unsigned int i=0; i<n_q_points; ++i)
-// psi_square[i] = sqr_point(psi[i]);
-
-// // add seminorm to L_2 norm or
-// // to zero
-// diff += inner_product (psi_square.begin(), psi_square.end(),
-// fe_values.get_JxW_values().begin(),
-// 0.0);
-// diff = sqrt(diff);
-
-// break;
-// };
+ // loop over all cells
+ DoFHandler<dim>::active_cell_iterator cell = dof.begin_active(),
+ endc = dof.end();
+ for (unsigned int index=0; cell != endc; ++cell, ++index)
+ {
+ double diff=0;
+ // initialize for this cell
+ fe_values.reinit (cell, boundary);
+
+ switch (norm)
+ {
+ case mean:
+ case L1_norm:
+ case L2_norm:
+ case Linfty_norm:
+ case H1_norm:
+ {
+ // we need the finite element
+ // function \psi at the different
+ // integration points. Compute
+ // it like this:
+ // \psi(x_j)=\sum_i v_i \phi_i(x_j)
+ // with v_i the nodal values of the
+ // fe_function and \phi_i(x_j) the
+ // matrix of the trial function
+ // values at the integration point
+ // x_j. Then the vector
+ // of the \psi(x_j) is v*Phi with
+ // v being the vector of nodal
+ // values on this cell and Phi
+ // the matrix.
+ //
+ // we then need the difference:
+ // reference_function(x_j)-\psi_j
+ // and assign that to the vector
+ // \psi.
+ const unsigned int n_q_points = q.n_quadrature_points;
+ vector<Vector<double> > psi (n_q_points);
+
+ // in praxi: first compute
+ // exact fe_function vector
+ exact_solution.value_list (fe_values.get_quadrature_points(),
+ psi);
+ // then subtract finite element
+ // fe_function
+ if (true)
+ {
+ vector< Vector<double> > function_values (n_q_points,
+ Vector<double>(fe.n_components));
+ fe_values.get_function_values (fe_function, function_values);
+
+/* transform (psi.begin(), psi.end(),
+ function_values.begin(),
+ psi.begin(),
+ minus<double>());
+*/ };
+
+ // for L1_norm and Linfty_norm:
+ // take absolute
+ // value, for the L2_norm take
+ // square of psi
+/* switch (norm)
+ {
+ case mean:
+ break;
+ case L1_norm:
+ case Linfty_norm:
+ transform (psi.begin(), psi.end(),
+ psi.begin(), ptr_fun(fabs));
+ break;
+ case L2_norm:
+ case H1_norm:
+ transform (psi.begin(), psi.end(),
+ psi.begin(), ptr_fun(sqr));
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+*/
+ // ok, now we have the integrand,
+ // let's compute the integral,
+ // which is
+ // sum_j psi_j JxW_j
+ // (or |psi_j| or |psi_j|^2
+/* switch (norm)
+ {
+ case mean:
+ case L1_norm:
+ diff = inner_product (psi.begin(), psi.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0);
+ break;
+ case L2_norm:
+ case H1_norm:
+ diff = sqrt(inner_product (psi.begin(), psi.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0));
+ break;
+ case Linfty_norm:
+ diff = *max_element (psi.begin(), psi.end());
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+*/
+ // note: the H1_norm uses the result
+ // of the L2_norm and control goes
+ // over to the next case statement!
+ if (norm != H1_norm)
+ break;
+ };
+
+ case H1_seminorm:
+ {
+ // note: the computation of the
+ // H1_norm starts at the previous
+ // case statement, but continues
+ // here!
+
+ // for H1_norm: re-square L2_norm.
+ diff = sqr(diff);
+
+ // same procedure as above, but now
+ // psi is a vector of gradients
+ const unsigned int n_q_points = q.n_quadrature_points;
+ vector<Tensor<1,dim> > psi (n_q_points);
+
+ // in praxi: first compute
+ // exact fe_function vector
+/* exact_solution.gradient_list (fe_values.get_quadrature_points(),
+ psi);
+*/
+ // then subtract finite element
+ // fe_function
+ if (true)
+ {
+ vector<Tensor<1,dim> > function_grads (n_q_points, Tensor<1,dim>());
+ fe_values.get_function_grads (fe_function, function_grads);
+
+/* transform (psi.begin(), psi.end(),
+ function_grads.begin(),
+ psi.begin(),
+ minus<Tensor<1,dim> >());
+*/ };
+ // take square of integrand
+ vector<double> psi_square (psi.size(), 0.0);
+ for (unsigned int i=0; i<n_q_points; ++i)
+ psi_square[i] = sqr_point(psi[i]);
+
+ // add seminorm to L_2 norm or
+ // to zero
+/* diff += inner_product (psi_square.begin(), psi_square.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0);
+ diff = sqrt(diff);
+*/
+ break;
+ };
-// default:
-// Assert (false, ExcNotImplemented());
-// };
+ default:
+ Assert (false, ExcNotImplemented());
+ };
-// // append result of this cell
-// // to the end of the vector
-// difference(index) = diff;
-// };
-// };
+ // append result of this cell
+ // to the end of the vector
+ difference(index) = diff;
+ };
+ };
template VectorTools<deal_II_dimension>;
{
next if (m/DeclException/);
chop;
- if (m/(.+)=(.+)/)
+ s/\</\<\;/g;
+ s/\>/\>\;/g;
+ if (m/(.+)=([^=]+)=?$/)
{
@n = split "::", $1;
{
$entry .= "+" . pop @n;
}
- $entry .= '=' . $2 . '=' . $library;
+ $entry .= ',' . $2 . ',' . $library;
push @entries, $entry;
}
else
;
foreach $entry (sort @entries)
{
- @l = split "=", $entry;
+ @l = split ",", $entry;
@c = split '\+', @l[0];
@f = split "#", @l[1];