typename DH::active_cell_iterator cell = cell_hint.get();
if (cell == dh->end())
cell = dh->begin_active();
- Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
- // Check if we already have all we need
- if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+ boost::optional<Point<dim> >
+ qp = get_reference_coordinates (cell, p);
+ if (!qp)
{
const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
= GridTools::find_active_cell_around_point (mapping, *dh, p);
cell_hint.get() = cell;
// Now we can find out about the point
- Quadrature<dim> quad(qp);
+ Quadrature<dim> quad(qp.get());
FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
update_values);
fe_v.reinit(cell);
}
+
template <int dim, typename DH, typename VECTOR>
void
- FEFieldFunction<dim, DH, VECTOR>::vector_gradient
- (const Point<dim> &p,
- std::vector<Tensor<1,dim> > &gradients) const
+ FEFieldFunction<dim, DH, VECTOR>::
+ vector_gradient (const Point<dim> &p,
+ std::vector<Tensor<1,dim> > &gradients) const
{
Assert (gradients.size() == n_components,
ExcDimensionMismatch(gradients.size(), n_components));
typename DH::active_cell_iterator cell = cell_hint.get();
if(cell == dh->end())
cell = dh->begin_active();
- Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
- // Check if we already have all we need
- if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+ boost::optional<Point<dim> >
+ qp = get_reference_coordinates (cell, p);
+ if (!qp)
{
std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
= GridTools::find_active_cell_around_point (mapping, *dh, p);
cell_hint.get() = cell;
// Now we can find out about the point
- Quadrature<dim> quad(qp);
+ Quadrature<dim> quad(qp.get());
FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
update_gradients);
fe_v.reinit(cell);
template <int dim, typename DH, typename VECTOR>
- Tensor<1,dim> FEFieldFunction<dim, DH, VECTOR>::gradient
- (const Point<dim> &p, unsigned int comp) const
+ Tensor<1,dim>
+ FEFieldFunction<dim, DH, VECTOR>::
+ gradient (const Point<dim> &p,
+ const unsigned int comp) const
{
std::vector<Tensor<1,dim> > grads(n_components);
vector_gradient(p, grads);
template <int dim, typename DH, typename VECTOR>
- void FEFieldFunction<dim, DH, VECTOR>::vector_laplacian
- (const Point<dim> &p,Vector<double> &values) const
+ void
+ FEFieldFunction<dim, DH, VECTOR>::
+ vector_laplacian (const Point<dim> &p,
+ Vector<double> &values) const
{
Assert (values.size() == n_components,
ExcDimensionMismatch(values.size(), n_components));
typename DH::active_cell_iterator cell = cell_hint.get();
if (cell == dh->end())
cell = dh->begin_active();
- Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
- // Check if we already have all we need
- if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+ boost::optional<Point<dim> >
+ qp = get_reference_coordinates (cell, p);
+ if (!qp)
{
const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
= GridTools::find_active_cell_around_point (mapping, *dh, p);
cell_hint.get() = cell;
// Now we can find out about the point
- Quadrature<dim> quad(qp);
+ Quadrature<dim> quad(qp.get());
FEValues<dim> fe_v(mapping, cell->get_fe(), quad,
update_hessians);
fe_v.reinit(cell);
points.push_back (Point<dim>(-0.7+i*0.07,-0.7+j*0.07));
points.push_back (Point<dim>(-0.27999999999999992, -0.62999999999999989));
- std::vector<Vector<double> > m (points.size());
+ std::vector<Vector<double> > m (points.size(), Vector<double>(2));
fe_function.vector_value_list (points, m);
for (unsigned int i=0; i<m.size(); ++i)
4,
&dof_handler));
- std::vector<Vector<double> > m (points.size());
+ std::vector<Vector<double> > m (points.size(), Vector<double>(2));
fe_function.vector_value_list (points, m);
for (unsigned int i=0; i<m.size(); ++i)
--- /dev/null
+//-------------------------------------------------------
+// $Id: fe_field_function_05_vector.cc $
+// Version: $Name$
+//
+// Copyright (C) 2005, 2011, 2012 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.
+//
+//-------------------------------------------------------
+
+
+// FEFieldFunction ran into an assertion after
+// Mapping::transform_real_to_unit_cell started throwing exceptions
+// when it couldn't find the point on the reference cell that belongs
+// to a given point, rather than just silently giving up
+//
+// this is yet another variant of _04 and _05 found by inspecting the
+// code in FEFieldFunction but the (same) exception is triggered in a
+// different place
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vectors.h>
+#include <deal.II/fe/fe_system.h>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+ public:
+ F() : Function<dim>(2) {}
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ v = 0;
+ v[0] = p.square();
+ }
+};
+
+
+
+template<int dim>
+void test()
+{
+ const HyperBallBoundary<dim> boundary_description;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_ball (triangulation);
+ triangulation.set_boundary (0, boundary_description);
+ triangulation.refine_global (1);
+
+ FESystem<dim> fe(FE_Q<dim> (2),2);
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(fe);
+
+ // interpolate a quadratic function
+ // into the space; this function
+ // can be represented exactly, so
+ // that we can compare again later
+ Vector<double> solution(dof_handler.n_dofs());
+ VectorTools::interpolate (dof_handler, F<dim>(), solution);
+
+ Functions::FEFieldFunction<2> fe_function (dof_handler, solution);
+
+ // add only one points but also set
+ // the active cell to one that
+ // doesn't contain the current
+ // point. the problem happens
+ // because we walk over a bunch of
+ // cells in the process of finding
+ // all of these points and then
+ // realize when we get to the one
+ // at the end that the coordinates
+ // for this point can't be found in
+ // the cell we have touched last
+ // (it's too far away from that
+ // cell, and the inverse mapping
+ // does not converge
+ Point<dim> point(-0.27999999999999992, -0.62999999999999989);
+ fe_function.set_active_cell (typename DoFHandler<dim>::active_cell_iterator
+ (&triangulation,
+ 1,
+ 4,
+ &dof_handler));
+
+ Vector<double> m(2);
+ fe_function.vector_value (point, m);
+
+ {
+ Assert (std::fabs(m(0) - point.square())
+ <
+ 1e-10 * std::fabs(m(0) + point.square()),
+ ExcInternalError());
+
+ Assert (std::fabs(m(1))
+ <
+ 1e-10,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+ std::ofstream logfile("fe_field_function_06_vector/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test<2>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::OK
--- /dev/null
+//-------------------------------------------------------
+// $Id: fe_field_function_05_vector.cc $
+// Version: $Name$
+//
+// Copyright (C) 2005, 2011, 2012 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.
+//
+//-------------------------------------------------------
+
+
+// FEFieldFunction ran into an assertion after
+// Mapping::transform_real_to_unit_cell started throwing exceptions
+// when it couldn't find the point on the reference cell that belongs
+// to a given point, rather than just silently giving up
+//
+// this is yet another variant of _04 and _05 and _06 found by
+// inspecting the code in FEFieldFunction but the (same) exception is
+// triggered in a different place
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vectors.h>
+#include <deal.II/fe/fe_system.h>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+ public:
+ F() : Function<dim>(2) {}
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ v[0] = p[0];
+ v[1] = 0;
+ }
+
+ virtual void vector_gradient (const Point<dim> &p,
+ std::vector<Tensor<1,dim> > &v) const
+ {
+ v[0] = 0;
+ v[1] = 0;
+ v[0][0] = 1;
+ }
+};
+
+
+
+template<int dim>
+void test()
+{
+ const HyperBallBoundary<dim> boundary_description;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_ball (triangulation);
+ triangulation.set_boundary (0, boundary_description);
+ triangulation.refine_global (1);
+
+ FESystem<dim> fe(FE_Q<dim> (2),2);
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(fe);
+
+ // interpolate a quadratic function
+ // into the space; this function
+ // can be represented exactly, so
+ // that we can compare again later
+ Vector<double> solution(dof_handler.n_dofs());
+ VectorTools::interpolate (dof_handler, F<dim>(), solution);
+
+ Functions::FEFieldFunction<2> fe_function (dof_handler, solution);
+
+ // add only one points but also set
+ // the active cell to one that
+ // doesn't contain the current
+ // point. the problem happens
+ // because we walk over a bunch of
+ // cells in the process of finding
+ // all of these points and then
+ // realize when we get to the one
+ // at the end that the coordinates
+ // for this point can't be found in
+ // the cell we have touched last
+ // (it's too far away from that
+ // cell, and the inverse mapping
+ // does not converge
+ Point<dim> point(-0.27999999999999992, -0.62999999999999989);
+ fe_function.set_active_cell (typename DoFHandler<dim>::active_cell_iterator
+ (&triangulation,
+ 1,
+ 4,
+ &dof_handler));
+
+ std::vector<Tensor<1,dim> > m(2);
+ fe_function.vector_gradient (point, m);
+
+ Assert (std::fabs(m[0][0] - 1)
+ <
+ 1e-10 * std::fabs(m[0][0] + 1),
+ ExcInternalError());
+ Assert (std::fabs(m[0][1])
+ <
+ 1e-10,
+ ExcInternalError());
+ Assert (m[1].norm()
+ <
+ 1e-10,
+ ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+ std::ofstream logfile("fe_field_function_07_vector/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test<2>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::OK
--- /dev/null
+//-------------------------------------------------------
+// $Id: fe_field_function_08_vector.cc $
+// Version: $Name$
+//
+// Copyright (C) 2005, 2011, 2012 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.
+//
+//-------------------------------------------------------
+
+
+// FEFieldFunction ran into an assertion after
+// Mapping::transform_real_to_unit_cell started throwing exceptions
+// when it couldn't find the point on the reference cell that belongs
+// to a given point, rather than just silently giving up
+//
+// this is yet another variant of _04 through _07 found by inspecting the
+// code in FEFieldFunction but the (same) exception is triggered in a
+// different place
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_q1.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/numerics/vectors.h>
+#include <deal.II/fe/fe_system.h>
+
+
+template <int dim>
+class F : public Function<dim>
+{
+ public:
+ F() : Function<dim>(2) {}
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ v = 0;
+ for (unsigned int i=0; i<dim; ++i)
+ v[0] += p[i]*p[i]*p[i]*p[i];
+ }
+
+ virtual void vector_laplacian (const Point<dim> &p,
+ Vector<double> &v) const
+ {
+ v = 0;
+ v[0] = p.square() * 4 * 3 * dim;
+ }
+};
+
+
+
+template<int dim>
+void test()
+{
+ const HyperBallBoundary<dim> boundary_description;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_ball (triangulation);
+ triangulation.set_boundary (0, boundary_description);
+ triangulation.refine_global (1);
+
+ FESystem<dim> fe(FE_Q<dim> (4),2);
+ DoFHandler<dim> dof_handler(triangulation);
+ dof_handler.distribute_dofs(fe);
+
+ // interpolate a quadratic function
+ // into the space; this function
+ // can be represented exactly, so
+ // that we can compare again later
+ Vector<double> solution(dof_handler.n_dofs());
+ VectorTools::interpolate (dof_handler, F<dim>(), solution);
+
+ Functions::FEFieldFunction<2> fe_function (dof_handler, solution);
+
+ // add only one points but also set
+ // the active cell to one that
+ // doesn't contain the current
+ // point. the problem happens
+ // because we walk over a bunch of
+ // cells in the process of finding
+ // all of these points and then
+ // realize when we get to the one
+ // at the end that the coordinates
+ // for this point can't be found in
+ // the cell we have touched last
+ // (it's too far away from that
+ // cell, and the inverse mapping
+ // does not converge
+ Point<dim> point(-0.27999999999999992, -0.62999999999999989);
+ fe_function.set_active_cell (typename DoFHandler<dim>::active_cell_iterator
+ (&triangulation,
+ 1,
+ 4,
+ &dof_handler));
+
+ Vector<double> m(2);
+ fe_function.vector_laplacian (point, m);
+
+ {
+ Assert (std::fabs(m(0) - point.square()*4*3)
+ <
+ 1e-8 * std::fabs(m(0) + point.square()*4*3),
+ ExcInternalError());
+
+ Assert (std::fabs(m(1))
+ <
+ 1e-10,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+ std::ofstream logfile("fe_field_function_08_vector/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test<2>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::OK
projects / dealii-svn.git / commitdiff