From 2e253b6fb4beacad76cf3d152ffa301b0def475c Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 26 Feb 1999 13:56:52 +0000 Subject: [PATCH] Errors according to Tensors and new LAC removed git-svn-id: https://svn.dealii.org/trunk@916 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/tensor.h | 28 +- deal.II/base/include/base/tensor_base.h | 17 +- deal.II/base/include/base/tensor_function.h | 14 +- deal.II/base/source/tensor.cc | 31 +- deal.II/base/source/tensor_function.cc | 7 +- deal.II/deal.II/include/fe/fe_values.h | 2 +- deal.II/deal.II/include/grid/tria.h | 12 +- deal.II/deal.II/include/numerics/vectors.h | 16 +- deal.II/deal.II/source/fe/fe_system.cc | 8 +- deal.II/deal.II/source/fe/fe_values.cc | 4 +- deal.II/deal.II/source/grid/tria.cc | 65 +++- deal.II/deal.II/source/numerics/data_io.cc | 13 +- deal.II/deal.II/source/numerics/vectors.cc | 366 ++++++++++---------- deal.II/doc/auto/scripts/index.pl | 8 +- 14 files changed, 329 insertions(+), 262 deletions(-) diff --git a/deal.II/base/include/base/tensor.h b/deal.II/base/include/base/tensor.h index a9ca5b8d21..93ae01a757 100644 --- a/deal.II/base/include/base/tensor.h +++ b/deal.II/base/include/base/tensor.h @@ -138,7 +138,7 @@ class Tensor { * usual in C++, the rightmost * index marches fastest. */ - void unroll(vector & result) const; + void unroll(Vector & result) const; /** @@ -146,18 +146,18 @@ class Tensor { */ 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 @@ -168,10 +168,10 @@ class Tensor { /** * Help function for unroll. */ - void unroll_recursion(vector & result, unsigned& start_index) const; + void unroll_recursion(Vector & result, unsigned& start_index) const; template<> - friend void Tensor::unroll_recursion(vector &, unsigned&) const; + friend void Tensor::unroll_recursion(Vector &, unsigned&) const; }; diff --git a/deal.II/base/include/base/tensor_base.h b/deal.II/base/include/base/tensor_base.h index af991fbed6..49b6fec21e 100644 --- a/deal.II/base/include/base/tensor_base.h +++ b/deal.II/base/include/base/tensor_base.h @@ -13,8 +13,7 @@ #include #include #include - - +#include // general template; specialized for rank==1; the general template is in @@ -171,7 +170,7 @@ class Tensor<1,dim> { * usual in C++, the rightmost * index marches fastest. */ - void unroll(vector & result) const; + void unroll(Vector & result) const; protected: /** @@ -182,15 +181,17 @@ class Tensor<1,dim> { /** * Help function for unroll. */ - void unroll_recursion(vector & result, unsigned& start_index) const; + void unroll_recursion(Vector & result, unsigned& start_index) const; template<> - friend void Tensor<2,dim>::unroll_recursion(vector &, unsigned&) const; + friend void Tensor<2,dim>::unroll_recursion(Vector &, 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); diff --git a/deal.II/base/include/base/tensor_function.h b/deal.II/base/include/base/tensor_function.h index f79825e2f2..810a2d6300 100644 --- a/deal.II/base/include/base/tensor_function.h +++ b/deal.II/base/include/base/tensor_function.h @@ -11,7 +11,7 @@ #include #include #include - +#include template class Tensor; @@ -53,7 +53,7 @@ class VectorFunction : * the same size as the #n_components# * array. */ - virtual void value (const Point &p, vector &values) const; + virtual void value (const Point &p, Vector &values) const; /** * Set #values# to the point values @@ -64,7 +64,7 @@ class VectorFunction : * array. */ virtual void value_list (const vector > &points, - vector > &values) const; + vector > &values) const; /** * Set #gradients# to the gradients of @@ -172,11 +172,17 @@ class TensorFunction : virtual void gradient_list (const vector > &points, vector > &gradients) const; + /** + * See #VectorFunction#. + */ + virtual void value (const Point &points, + Vector &values) const; + /** * See #VectorFunction#. */ virtual void value_list (const vector > &points, - vector > &values) const; + vector > &values) const; /** * See #VectorFunction#. diff --git a/deal.II/base/source/tensor.cc b/deal.II/base/source/tensor.cc index 3a4bac21a1..d1035b5008 100644 --- a/deal.II/base/source/tensor.cc +++ b/deal.II/base/source/tensor.cc @@ -2,17 +2,24 @@ #include #include +#include -template void -Tensor<1,dim>::unroll( vector& result) const +template +void +Tensor<1,dim>::unroll( Vector& result) const { - Assert(false, ExcNotImplemented()); + Assert(result.size()==dim, + ExcWrongVectorSize(dim, result.size())); + + unsigned index = 0; + unroll_recursion(result,index); } -template void -Tensor::unroll( vector& result) const +template +void +Tensor::unroll( Vector& result) const { Assert(result.size()==pow(dim,rank_), ExcWrongVectorSize(pow(dim,rank_), result.size())); @@ -22,8 +29,9 @@ Tensor::unroll( vector& result) const } -template void -Tensor::unroll_recursion( vector& result, unsigned& index) const +template +void +Tensor::unroll_recursion( Vector& result, unsigned& index) const { for (unsigned i=0; i::unroll_recursion( vector& result, unsigned& index) c } -template void -Tensor<1,dim>::unroll_recursion( vector& result, unsigned& index) const +template +void +Tensor<1,dim>::unroll_recursion( Vector& result, unsigned& index) const { for (unsigned i=0; i #include #include +#include template VectorFunction::VectorFunction(unsigned n_components, const double initial_time) @@ -29,7 +30,7 @@ VectorFunction::operator () (const Point &, unsigned) const */ template void -VectorFunction::value (const Point &, vector &) const +VectorFunction::value (const Point &, Vector &) const { Assert (false, ExcPureFunctionCalled()); } @@ -37,7 +38,7 @@ VectorFunction::value (const Point &, vector &) const template void VectorFunction::value_list (const vector > &, - vector > &) const + vector > &) const { Assert (false, ExcPureFunctionCalled()); } @@ -139,7 +140,7 @@ TensorFunction::value (const Point &p, vector &erg) co template void TensorFunction::value_list (const vector > & points, - vector > & values) const + vector > & values) const { Assert (values.size() == points.size(), ExcVectorHasWrongSize(values.size(), points.size())); diff --git a/deal.II/deal.II/include/fe/fe_values.h b/deal.II/deal.II/include/fe/fe_values.h index b57150d59f..d55dc026c2 100644 --- a/deal.II/deal.II/include/fe/fe_values.h +++ b/deal.II/deal.II/include/fe/fe_values.h @@ -304,7 +304,7 @@ class FEValuesBase { * elements. */ void get_function_values (const Vector &fe_function, - vector > &values) const; + vector > &values) const; /** * Return the gradient of the #i#th shape diff --git a/deal.II/deal.II/include/grid/tria.h b/deal.II/deal.II/include/grid/tria.h index e3d6de9137..cddb2210c2 100644 --- a/deal.II/deal.II/include/grid/tria.h +++ b/deal.II/deal.II/include/grid/tria.h @@ -1676,7 +1676,8 @@ class Triangulation : public TriaDimensionInfo, public Subscriptor { * much needed here and saving memory may * be a good goal when using many cells. */ - void refine (const Vector &criteria, + template + void refine (const Vector &criteria, const double threshold); /** @@ -1692,7 +1693,8 @@ class Triangulation : public TriaDimensionInfo, public Subscriptor { * much needed here and saving memory may * be a good goal when using many cells. */ - void coarsen (const Vector &criteria, + template + void coarsen (const Vector &criteria, const double threshold); /** @@ -1716,7 +1718,8 @@ class Triangulation : public TriaDimensionInfo, public Subscriptor { * much needed here and saving memory may * be a good goal when using many cells. */ - void refine_and_coarsen_fixed_number (const Vector &criteria, + template + void refine_and_coarsen_fixed_number (const Vector &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells); @@ -1741,7 +1744,8 @@ class Triangulation : public TriaDimensionInfo, public Subscriptor { * much needed here and saving memory may * be a good goal when using many cells. */ - void refine_and_coarsen_fixed_fraction (const Vector &criteria, + template + void refine_and_coarsen_fixed_fraction (const Vector &criteria, const double top_fraction, const double bottom_fraction); diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h index a4567f8a75..f3b9c6e57c 100644 --- a/deal.II/deal.II/include/numerics/vectors.h +++ b/deal.II/deal.II/include/numerics/vectors.h @@ -398,14 +398,14 @@ class VectorTools { * Compute the error for the solution of a system. * See the other #integrate_difference#. */ -// static void integrate_difference (const DoFHandler &dof, -// const Vector &fe_function, -// const TensorFunction<1,dim>&exact_solution, -// Vector &difference, -// const Quadrature &q, -// const FiniteElement &fe, -// const NormType &norm, -// const Boundary &boundary); + static void integrate_difference (const DoFHandler &dof, + const Vector &fe_function, + const VectorFunction&exact_solution, + Vector &difference, + const Quadrature &q, + const FiniteElement &fe, + const NormType &norm, + const Boundary &boundary); /** diff --git a/deal.II/deal.II/source/fe/fe_system.cc b/deal.II/deal.II/source/fe/fe_system.cc index b9ae5a99b2..4415ea60f7 100644 --- a/deal.II/deal.II/source/fe/fe_system.cc +++ b/deal.II/deal.II/source/fe/fe_system.cc @@ -36,7 +36,7 @@ void FESystem::initialize () // 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::vertices_per_cell ; ++vertex_number) { unsigned comp_start = 0; @@ -60,7 +60,7 @@ void FESystem::initialize () } // 2. Lines - for (unsigned line_number= 0 ; line_number < GeometryInfo<2>::lines_per_cell ; + for (unsigned line_number= 0 ; line_number < GeometryInfo::lines_per_cell ; ++line_number) { unsigned comp_start = 0; @@ -85,7 +85,7 @@ void FESystem::initialize () } // 3. Quads - for (unsigned quad_number= 0 ; quad_number < GeometryInfo<2>::quads_per_cell ; + for (unsigned quad_number= 0 ; quad_number < GeometryInfo::quads_per_cell ; ++quad_number) { unsigned comp_start = 0; @@ -110,7 +110,7 @@ void FESystem::initialize () } // 4. Hex - for (unsigned hex_number= 0 ; hex_number < GeometryInfo<2>::hexes_per_cell ; + for (unsigned hex_number= 0 ; hex_number < GeometryInfo::hexes_per_cell ; ++hex_number) { unsigned comp_start = 0; diff --git a/deal.II/deal.II/source/fe/fe_values.cc b/deal.II/deal.II/source/fe/fe_values.cc index 67c57e16fa..31568f053c 100644 --- a/deal.II/deal.II/source/fe/fe_values.cc +++ b/deal.II/deal.II/source/fe/fe_values.cc @@ -90,7 +90,7 @@ void FEValuesBase::get_function_values (const Vector &fe_function, template void FEValuesBase::get_function_values (const Vector &fe_function, - vector< vector > &values) const + vector< Vector > &values) const { Assert (fe->n_components == values.size(), ExcWrongNoOfComponents()); @@ -113,7 +113,7 @@ void FEValuesBase::get_function_values (const Vector &fe_functi // functions for (unsigned int point=0; pointsystem_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)); }; diff --git a/deal.II/deal.II/source/grid/tria.cc b/deal.II/deal.II/source/grid/tria.cc index 605ac44c11..85108aead8 100644 --- a/deal.II/deal.II/source/grid/tria.cc +++ b/deal.II/deal.II/source/grid/tria.cc @@ -3761,9 +3761,10 @@ void Triangulation<3>::print_gnuplot (ostream &out, -template -void Triangulation::refine (const Vector &criteria, - const double threshold) { +template template +void Triangulation::refine (const Vector &criteria, + const double threshold) +{ Assert (criteria.size() == n_active_cells(), ExcInvalidVectorSize(criteria.size(), n_active_cells())); @@ -3785,9 +3786,10 @@ void Triangulation::refine (const Vector &criteria, -template -void Triangulation::coarsen (const Vector &criteria, - const double threshold) { +template template +void Triangulation::coarsen (const Vector &criteria, + const double threshold) +{ Assert (criteria.size() == n_active_cells(), ExcInvalidVectorSize(criteria.size(), n_active_cells())); @@ -3801,10 +3803,12 @@ void Triangulation::coarsen (const Vector &criteria, -template -void Triangulation::refine_and_coarsen_fixed_number (const Vector &criteria, - const double top_fraction, - const double bottom_fraction) { +template template +void +Triangulation::refine_and_coarsen_fixed_number (const Vector &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()); @@ -3817,7 +3821,7 @@ void Triangulation::refine_and_coarsen_fixed_number (const Vector &c const int coarsen_cells = max(static_cast(bottom_fraction*criteria.size()), 1); - Vector tmp(criteria); + Vector tmp(criteria); nth_element (tmp.begin(), tmp.begin()+refine_cells, tmp.end(), greater()); @@ -3839,9 +3843,9 @@ double sqr(double a) { -template +template template void -Triangulation::refine_and_coarsen_fixed_fraction (const Vector &criteria, +Triangulation::refine_and_coarsen_fixed_fraction (const Vector &criteria, const double top_fraction, const double bottom_fraction) { // correct number of cells is @@ -3854,7 +3858,7 @@ Triangulation::refine_and_coarsen_fixed_fraction (const Vector &crit // error, which is what we have to sum // up and compare with // #fraction_of_error*total_error#. - Vector tmp(criteria); + Vector tmp(criteria); const double total_error = tmp.l1_norm(); Vector partial_sums(criteria.size()); @@ -7048,3 +7052,36 @@ void Triangulation::read_bool_vector (const unsigned int magic_number1, // explicit instantiations template class Triangulation; +template void Triangulation +::refine (const Vector &, const double); + +template void Triangulation +::refine (const Vector &, const double); + +template void Triangulation +::coarsen (const Vector &, const double); + +template void Triangulation +::coarsen (const Vector &, const double); + + +template void Triangulation +::refine_and_coarsen_fixed_number (const Vector &, + const double top_fraction, + const double bottom_fraction); + +template void Triangulation +::refine_and_coarsen_fixed_number (const Vector &criteria, + const double top_fraction, + const double bottom_fraction); + +template void Triangulation +::refine_and_coarsen_fixed_fraction (const Vector &criteria, + const double top_fraction, + const double bottom_fraction); + +template void Triangulation +::refine_and_coarsen_fixed_fraction (const Vector &criteria, + const double top_fraction, + const double bottom_fraction); + diff --git a/deal.II/deal.II/source/numerics/data_io.cc b/deal.II/deal.II/source/numerics/data_io.cc index da9a132577..b34ce01796 100644 --- a/deal.II/deal.II/source/numerics/data_io.cc +++ b/deal.II/deal.II/source/numerics/data_io.cc @@ -543,8 +543,11 @@ void DataOut::write_gnuplot (ostream &out, unsigned int accuracy) const FEValues fe(dofs->get_fe(), points, UpdateFlags(update_q_points)); const StraightBoundary boundary; - vector< vector > > - values (data.size(), vector< vector >(dofs->get_fe().n_components, vector(points.n_quadrature_points))); + vector< vector > > + values (data.size(), + vector< Vector >(points.n_quadrature_points, + Vector(dofs->get_fe().n_components + ))); for (cell=dofs->begin_active(); cell!=endc; ++cell) { @@ -570,7 +573,7 @@ void DataOut::write_gnuplot (ostream &out, unsigned int accuracy) const 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; }; @@ -592,7 +595,7 @@ void DataOut::write_gnuplot (ostream &out, unsigned int accuracy) const 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; } @@ -625,7 +628,7 @@ void DataOut::write_gnuplot (ostream &out, unsigned int accuracy) const 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; } diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc index 31330e7b62..68785d54dd 100644 --- a/deal.II/deal.II/source/numerics/vectors.cc +++ b/deal.II/deal.II/source/numerics/vectors.cc @@ -4,6 +4,7 @@ #include +#include #include #include #include @@ -557,193 +558,196 @@ void VectorTools::integrate_difference (const DoFHandler &dof, -// template -// void VectorTools::integrate_difference (const DoFHandler &dof, -// const Vector &fe_function, -// const TensorFunction<1,dim>&exact_solution, -// Vector &difference, -// const Quadrature &q, -// const FiniteElement &fe, -// const NormType &norm, -// const Boundary &boundary) { -// Assert (fe == dof.get_fe(), ExcInvalidFE()); +template +void +VectorTools::integrate_difference (const DoFHandler &dof, + const Vector &fe_function, + const VectorFunction&exact_solution, + Vector &difference, + const Quadrature &q, + const FiniteElement &fe, + const NormType &norm, + const Boundary &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 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 fe_values(fe, q, update_flags); -// // loop over all cells -// DoFHandler::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 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 > function_values (fe.n_components, vector(n_q_points, 0)); -// fe_values.get_function_values (fe_function, function_values); - -// transform (psi.begin(), psi.end(), -// function_values.begin(), -// psi.begin(), -// minus()); -// }; - -// // 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 > 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 > 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 >()); -// }; -// // take square of integrand -// vector psi_square (psi.size(), 0.0); -// for (unsigned int i=0; i::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 > 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 > function_values (n_q_points, + Vector(fe.n_components)); + fe_values.get_function_values (fe_function, function_values); + +/* transform (psi.begin(), psi.end(), + function_values.begin(), + psi.begin(), + minus()); +*/ }; + + // 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 > 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 > 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 >()); +*/ }; + // take square of integrand + vector psi_square (psi.size(), 0.0); + for (unsigned int i=0; i; diff --git a/deal.II/doc/auto/scripts/index.pl b/deal.II/doc/auto/scripts/index.pl index 737c091217..85a4086f55 100644 --- a/deal.II/doc/auto/scripts/index.pl +++ b/deal.II/doc/auto/scripts/index.pl @@ -4,7 +4,9 @@ while(<>) { next if (m/DeclException/); chop; - if (m/(.+)=(.+)/) + s/\/\>\;/g; + if (m/(.+)=([^=]+)=?$/) { @n = split "::", $1; @@ -13,7 +15,7 @@ while(<>) { $entry .= "+" . pop @n; } - $entry .= '=' . $2 . '=' . $library; + $entry .= ',' . $2 . ',' . $library; push @entries, $entry; } else @@ -34,7 +36,7 @@ EOT ; foreach $entry (sort @entries) { - @l = split "=", $entry; + @l = split ",", $entry; @c = split '\+', @l[0]; @f = split "#", @l[1]; -- 2.39.5