* simple test for checking
* whether a point is in a cell,
* it is only implemented for
- * Q1-mapping yet.
+ * Q1-mapping at this time.
*/
template <int dim, class InVector>
static void point_difference (const DoFHandler<dim>& dof,
Vector<double>& difference,
const Point<dim>& point);
+ /**
+ * Evaluate a possibly
+ * vector-valued finite element
+ * function defined by the given
+ * DoFHandler and nodal vector at
+ * the given point, and return
+ * the (vector) value of this
+ * function through the last
+ * argument.
+ */
+ template <int dim, class InVector>
+ static
+ void
+ point_value (const DoFHandler<dim> &dof,
+ const InVector &fe_function,
+ const Point<dim> &point,
+ Vector<double> &value);
+
+ /**
+ * Evaluate a scalar finite
+ * element function defined by
+ * the given DoFHandler and nodal
+ * vector at the given point, and
+ * return the value of this
+ * function.
+ */
+ template <int dim, class InVector>
+ static
+ double
+ point_value (const DoFHandler<dim> &dof,
+ const InVector &fe_function,
+ const Point<dim> &point);
+
/**
* Subtract the (algebraic) mean value
* from a vector. This function is most
// $Id$
// Version: $Name$
//
-// Copyright (C) 2005 by the deal.II authors
+// Copyright (C) 2005, 2006 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
ExcInternalError());
- const Quadrature<dim> quadrature (std::vector<Point<dim> > (1, unit_point),
- std::vector<double> (1, 1.));
+ const Quadrature<dim> quadrature (unit_point);
FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
fe_values.reinit(cell);
}
+template <int dim, class InVector>
+void
+VectorTools::point_value (const DoFHandler<dim> &dof,
+ const InVector &fe_function,
+ const Point<dim> &point,
+ Vector<double> &value)
+{
+ static const MappingQ1<dim> mapping;
+ const FiniteElement<dim>& fe = dof.get_fe();
+
+ Assert(value.size() == fe.n_components(),
+ ExcDimensionMismatch(value.size(), fe.n_components()));
+
+ // first find the cell in which this point
+ // is, initialize a quadrature rule with
+ // it, and then a FEValues object
+ const typename DoFHandler<dim>::active_cell_iterator
+ cell = GridTools::find_active_cell_around_point (dof, point);
+
+ const Point<dim> unit_point
+ = mapping.transform_real_to_unit_cell(cell, point);
+ Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
+ ExcInternalError());
+
+ const Quadrature<dim> quadrature (unit_point);
+ FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+ fe_values.reinit(cell);
+
+ // then use this to get at the values of
+ // the given fe_function at this point
+ std::vector<Vector<double> > u_value(1, Vector<double> (fe.n_components()));
+ fe_values.get_function_values(fe_function, u_value);
+
+ value = u_value[0];
+}
+
+
+
+template <int dim, class InVector>
+double
+VectorTools::point_value (const DoFHandler<dim> &dof,
+ const InVector &fe_function,
+ const Point<dim> &point)
+{
+ static const MappingQ1<dim> mapping;
+ const FiniteElement<dim>& fe = dof.get_fe();
+
+ Assert(fe.n_components() == 1,
+ ExcMessage ("Finite element is not scalar as is necessary for this function"));
+
+ // first find the cell in which this point
+ // is, initialize a quadrature rule with
+ // it, and then a FEValues object
+ const typename DoFHandler<dim>::active_cell_iterator
+ cell = GridTools::find_active_cell_around_point (dof, point);
+
+ const Point<dim> unit_point
+ = mapping.transform_real_to_unit_cell(cell, point);
+ Assert (GeometryInfo<dim>::is_inside_unit_cell (unit_point),
+ ExcInternalError());
+
+ const Quadrature<dim> quadrature (unit_point);
+ FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
+ fe_values.reinit(cell);
+
+ // then use this to get at the values of
+ // the given fe_function at this point
+ std::vector<double> u_value(1);
+ fe_values.get_function_values(fe_function, u_value);
+
+ return u_value[0];
+}
+
+
template <int dim, class InVector>
double
// $Id$
// Version: $Name$
//
-// Copyright (C) 2005 by the deal.II authors
+// Copyright (C) 2005, 2006 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
Vector<double>&,
const Point<deal_II_dimension>&);
+template
+void VectorTools::point_value<deal_II_dimension> (
+ const DoFHandler<deal_II_dimension>&,
+ const VEC&,
+ const Point<deal_II_dimension>&,
+ Vector<double>&);
+
+template
+double VectorTools::point_value<deal_II_dimension> (
+ const DoFHandler<deal_II_dimension>&,
+ const VEC&,
+ const Point<deal_II_dimension>&);
+
template
double VectorTools::compute_mean_value<deal_II_dimension>
(const Mapping<deal_II_dimension>&,
<h3>deal.II</h3>
<ol>
+ <li> <p>
+ New: There are now functions <code class="class">VectorTools</code>::<code
+ class="member">point_value</code> that evaluate the value of a
+ finite element function at a given point inside the domain.
+ <br>
+ (WB 2006/03/06)
+ </p>
<li> <p> Improved: <code class="class">GridOut</code> now has functions for
- declaring and parsing parameters in a <code
- class="class">ParameterHandler</code>, like used in <code
- class="class">DataOut_Interface</code> for a long time.
- <br>
- (GK 2006/03/06)
- </p>
+ declaring and parsing parameters in a <code
+ class="class">ParameterHandler</code>, like used in <code
+ class="class">DataOut_Interface</code> for a long time.
+ <br>
+ (GK 2006/03/06)
+ </p>
<li> <p>
Improved: The <code class="class">Triangulation</code>, <code
############################################################
# Makefile,v 1.15 2002/06/13 12:51:13 hartmann Exp
-# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors
+# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors
############################################################
############################################################
cylinder_shell_* \
compressed_sparsity_pattern_* \
point_difference_* \
+ point_value_* \
fe_q_3d* \
rt_* \
nedelec_* \
--- /dev/null
+//---------------------------- point_value_01.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- point_value_01.cc ---------------------------
+
+// check VectorTools::point_value
+
+
+
+#include "../tests.h"
+#include <base/logstream.h>
+#include <base/function_lib.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_accessor.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <fe/fe_q.h>
+#include <numerics/vectors.h>
+
+#include <fstream>
+#include <cmath>
+#include <iostream>
+
+
+template<int dim>
+class MySquareFunction : public Function<dim>
+{
+ public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const
+ { return (component+1)*p.square()+1; };
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+ { values(0) = value(p,0); };
+};
+
+
+template<int dim>
+class MyExpFunction : public Function<dim>
+{
+ public:
+ virtual double value (const Point<dim> &p,
+ const unsigned int component) const
+ { return std::exp (p(0)); };
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+ { values(0) = value(p,0); };
+};
+
+
+
+template <int dim>
+void make_mesh (Triangulation<dim> &tria)
+{
+
+ GridGenerator::hyper_cube(tria, -1, 1);
+
+ // refine the mesh in a random way so as to
+ // generate as many cells with
+ // hanging nodes as possible
+ tria.refine_global (4-dim);
+ const double steps[4] = { /*d=0*/ 0, 7, 3, 3 };
+ for (unsigned int i=0; i<steps[dim]; ++i)
+ {
+ typename Triangulation<dim>::active_cell_iterator
+ cell = tria.begin_active();
+ for (unsigned int index=0; cell != tria.end(); ++cell, ++index)
+ if (index % (3*dim) == 0)
+ cell->set_refine_flag();
+ tria.execute_coarsening_and_refinement ();
+ }
+}
+
+
+
+
+template <int dim>
+void
+check ()
+{
+ Triangulation<dim> tria;
+ make_mesh (tria);
+
+ FE_Q<dim> element(3);
+ DoFHandler<dim> dof(tria);
+ dof.distribute_dofs(element);
+
+ // test with two different functions: one
+ // that is exactly representable on the
+ // chosen finite element space, and one
+ // that isn't
+ for (unsigned int i=0; i<2; ++i)
+ {
+ static const MySquareFunction<dim> function_1;
+ static const Functions::CosineFunction<dim> function_2;
+
+ const Function<dim> &
+ function = (i==0 ?
+ static_cast<const Function<dim>&>(function_1) :
+ static_cast<const Function<dim>&>(function_2));
+
+ Vector<double> v (dof.n_dofs());
+ VectorTools::interpolate (dof, function, v);
+
+ // for the following points, check the
+ // function value, output it, and
+ // verify that the value retrieved from
+ // the interpolated function is close
+ // enough to that of the real function
+ //
+ // also verify that the actual value is
+ // roughly correct
+ Point<dim> p[3];
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ p[0][d] = 0;
+ p[1][d] = 0.5;
+ p[2][d] = 1./3.;
+ }
+ Vector<double> value(1);
+ for (unsigned int i=0; i<3; ++i)
+ {
+ VectorTools::point_value (dof, v, p[i], value);
+ deallog << -value(0) << std::endl;
+
+ Assert (std::abs(value(0) - function.value(p[i])) < 1e-4,
+ ExcInternalError());
+
+ const double scalar_value = VectorTools::point_value (dof, v, p[i]);
+ Assert (std::abs(value(0) - scalar_value) < 1e-4,
+ ExcInternalError());
+ }
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+ std::ofstream logfile ("point_value_01/output");
+ logfile.precision (4);
+ deallog.attach(logfile);
+ deallog.depth_console (0);
+
+ deallog.push ("1d");
+ check<1> ();
+ deallog.pop ();
+ deallog.push ("2d");
+ check<2> ();
+ deallog.pop ();
+ deallog.push ("3d");
+ check<3> ();
+ deallog.pop ();
+}
--- /dev/null
+
+DEAL:1d::-1.000
+DEAL:1d::-1.250
+DEAL:1d::-1.111
+DEAL:1d::-1.000
+DEAL:1d::-0.7071
+DEAL:1d::-0.8660
+DEAL:1d::OK
+DEAL:2d::-1.000
+DEAL:2d::-1.500
+DEAL:2d::-1.222
+DEAL:2d::-1.000
+DEAL:2d::-0.5000
+DEAL:2d::-0.7500
+DEAL:2d::OK
+DEAL:3d::-1.000
+DEAL:3d::-1.750
+DEAL:3d::-1.333
+DEAL:3d::-1.000
+DEAL:3d::-0.3536
+DEAL:3d::-0.6495
+DEAL:3d::OK