--- /dev/null
+//---------------------------- template.cc ---------------------------
+// $Id$
+// 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.
+//
+//---------------------------- template.cc ---------------------------
+
+
+// Test the FEFieldFunction class
+
+#include "../tests.h"
+#include <fstream>
+
+// all include files you need here
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vectors.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();
+ }
+};
+
+double abs_zero(double a)
+{
+ if( std::abs(a) < 1e-10)
+ return 0;
+ else
+ return a;
+}
+
+template <int dim>
+void test() {
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(8/dim);
+
+ FESystem<dim> fe(FE_Q<dim> (2),2);
+ DoFHandler<dim> dh(tria);
+
+ dh.distribute_dofs(fe);
+ deallog << "Ndofs :" << dh.n_dofs() << std::endl;
+
+
+
+ Vector<double> v1(dh.n_dofs()), v2(dh.n_dofs());
+
+ VectorTools::interpolate(dh, F<dim>(), v1);
+ deallog << "V norm: " << v1.l2_norm() << std::endl;
+
+ Functions::FEFieldFunction<dim, DoFHandler<dim>, Vector<double> >
+ fef(dh, v1);
+
+ VectorTools::interpolate(dh, fef, v2);
+
+ v2.add(-1, v1);
+ deallog << "Interpolation error: " << abs_zero(v2.l2_norm())
+ << std::endl;
+
+ Vector<double> error(tria.n_active_cells());
+ QGauss<dim> quad(2);
+ VectorTools::integrate_difference(dh, v1, F<dim>(),
+ error, quad,
+ VectorTools::L2_norm);
+ deallog << "L2 Interpolation error: "
+ << abs_zero(error.l2_norm()) << std::endl;
+ error = 0;
+ VectorTools::integrate_difference(dh, v1, fef,
+ error, quad,
+ VectorTools::L2_norm);
+ deallog << "L2 Interpolation error with fef: "
+ << abs_zero(error.l2_norm()) << std::endl;
+
+}
+
+int main ()
+{
+ std::ofstream logfile("fe_field_function_01_vector/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test<1>();
+ test<2>();
+ test<3>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Ndofs :1026
+DEAL::V norm: 10.1440
+DEAL::Interpolation error: 0.00000
+DEAL::L2 Interpolation error: 0.00000
+DEAL::L2 Interpolation error with fef: 0.00000
+DEAL::Ndofs :2178
+DEAL::V norm: 26.5652
+DEAL::Interpolation error: 0.00000
+DEAL::L2 Interpolation error: 0.00000
+DEAL::L2 Interpolation error with fef: 0.00000
+DEAL::Ndofs :1458
+DEAL::V norm: 32.6964
+DEAL::Interpolation error: 0.00000
+DEAL::L2 Interpolation error: 0.00000
+DEAL::L2 Interpolation error with fef: 0.00000
--- /dev/null
+//---------------------------- template.cc ---------------------------
+// $Id$
+// 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.
+//
+//---------------------------- template.cc ---------------------------
+
+
+// Test the FEFieldFunction class. This class was not thread-safe at one point
+// because it keeps a cache on the side that was invalidated when different
+// threads kept pouncing on it. Patrick Sodre wrote a fix for that that's
+// based on a thread-local variable instead of the single cache.
+
+#include "../tests.h"
+#include <fstream>
+
+// all include files you need here
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.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();
+ }
+};
+
+
+double abs_zero(double a)
+{
+ if( std::abs(a) < 1e-10)
+ return 0;
+ else
+ return a;
+}
+
+template <int dim>
+void test() {
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(8/dim);
+
+ FESystem<dim> fe(FE_Q<dim> (2),2);
+ DoFHandler<dim> dh(tria);
+
+ dh.distribute_dofs(fe);
+ deallog << "Ndofs :" << dh.n_dofs() << std::endl;
+
+ F<dim> ff;
+
+ Vector<double> v1(dh.n_dofs()), v2(dh.n_dofs());
+
+ VectorTools::interpolate(dh, ff, v1);
+ deallog << "V norm: " << v1.l2_norm() << std::endl;
+
+ Functions::FEFieldFunction<dim, DoFHandler<dim>, Vector<double> >
+ fef(dh, v1);
+
+ // project the discrete function fef back
+ // onto the finite element space. this
+ // should be the identity operation and
+ // consequently subtracting the vector from
+ // itself should yield zero
+ {
+ ConstraintMatrix cm;
+ cm.close();
+ VectorTools::project(dh, cm, QGauss<dim>(3), fef, v2);
+ }
+ v2.add(-1, v1);
+
+ deallog << "Projection error: " << abs_zero(v2.l2_norm())
+ << std::endl;
+ Assert (v2.l2_norm() < 1e-10, ExcInternalError());
+}
+
+int main ()
+{
+ std::ofstream logfile("fe_field_function_02_vector/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test<1>();
+ test<2>();
+ test<3>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Ndofs :1026
+DEAL::V norm: 10.1440
+DEAL::Projection error: 0.00000
+DEAL::Ndofs :2178
+DEAL::V norm: 26.5652
+DEAL::Projection error: 0.00000
+DEAL::Ndofs :1458
+DEAL::V norm: 32.6964
+DEAL::Projection error: 0.00000
--- /dev/null
+//-------------------------------------------------------
+// $Id: fe_field_function_03.cc $
+// Version: $Name$
+//
+// Copyright (C) 2005, 2011 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.
+//
+//-------------------------------------------------------
+
+
+// Test the functionality of the laplacian in the FEFieldFunction class.
+
+#include "../tests.h"
+#include <fstream>
+
+// all include files you need here
+#include <deal.II/numerics/fe_field_function.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/base/function_lib.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/numerics/vectors.h>
+
+template <int dim>
+void test() {
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(9/dim);
+
+ FE_Q<dim> fe(2);
+ DoFHandler<dim> dh(tria);
+
+ dh.distribute_dofs(fe);
+ deallog << "Ndofs :" << dh.n_dofs() << std::endl;
+
+ Functions::CosineFunction<dim> ff;
+
+ Vector<double> v1(dh.n_dofs()), v2(dh.n_dofs());
+
+ VectorTools::interpolate(dh, ff, v1);
+ deallog << "V norm: " << v1.l2_norm() << std::endl;
+
+ Functions::FEFieldFunction<dim, DoFHandler<dim>, Vector<double> >
+ fef(dh, v1);
+
+ //create the origin
+ Point<dim> p;
+ //compute the error of the laplacian in this point
+ deallog << "Value of the laplacian in 0:" << std::endl;
+ deallog << "correct value: " << ff.laplacian(p) <<", approximation: "<< fef.laplacian(p) << std::endl;
+
+ //now we want to test the list version
+ Point<dim> p1 = Point<dim>::unit_vector(0);
+ p1 = p1 * 0.5;
+ Point<dim> p2 = p1 * 0.5;
+ std::vector<Point<dim> > vec;
+ vec.push_back(p1);
+ vec.push_back(p2);
+ std::vector<double> values_c(2);
+ std::vector<double> values_a(2);
+
+ //get the laplacians at these two points
+ ff.laplacian_list(vec, values_c);
+ fef.laplacian_list(vec, values_a);
+ deallog << "Value of the laplacian in 0.5*e1 and 0.25 * e1:" << std::endl;
+ deallog << " correct values: " <<values_c[0] << " " << values_c[1] <<", approximations: "<<values_a[0] << " " << values_a[1] << std::endl;
+}
+
+int main ()
+{
+ std::ofstream logfile("fe_field_function_03/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test<1>();
+ test<2>();
+ test<3>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::Ndofs :1025
+DEAL::V norm: 22.6385
+DEAL::Value of the laplacian in 0:
+DEAL::correct value: -2.46740, approximation: -2.46743
+DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1:
+DEAL:: correct values: -1.74472 -2.27958, approximations: -1.74739 -2.28103
+DEAL::Ndofs :1089
+DEAL::V norm: 16.5000
+DEAL::Value of the laplacian in 0:
+DEAL::correct value: -4.93480, approximation: -4.92787
+DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1:
+DEAL:: correct values: -3.48943 -4.55916, approximations: -3.57012 -4.59908
+DEAL::Ndofs :4913
+DEAL::V norm: 24.7815
+DEAL::Value of the laplacian in 0:
+DEAL::correct value: -7.40220, approximation: -7.36064
+DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1:
+DEAL:: correct values: -5.23415 -6.83874, approximations: -5.37564 -6.89283
--- /dev/null
+//-------------------------------------------------------
+// $Id: fe_field_function_04_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
+
+#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);
+
+ Vector<double> solution(dof_handler.n_dofs());
+ VectorTools::interpolate (dof_handler, F<dim>(), solution);
+
+ Functions::FEFieldFunction<2> fe_function (dof_handler, solution);
+ std::vector<Point<dim> > points;
+
+ // add a bunch of points. all
+ // points are inside the circle.
+ // 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
+ for (unsigned int i=0; i<20; ++i)
+ for (unsigned int j=0; j<20; ++j)
+ 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());
+ fe_function.vector_value_list (points, m);
+
+ for (unsigned int i=0; i<m.size(); ++i)
+ {
+ Assert (std::fabs(m[i](0) - points[i].square())
+ <
+ 1e-10 * std::fabs(m[i](0) + points[i].square()),
+ ExcInternalError());
+
+ Assert (std::fabs(m[i](1))
+ <
+ 1e-10,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+ std::ofstream logfile("fe_field_function_04_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 a variant of _04 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);
+ std::vector<Point<dim> > points;
+
+ // 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
+ points.push_back (Point<dim>(-0.27999999999999992, -0.62999999999999989));
+ fe_function.set_active_cell (typename DoFHandler<dim>::active_cell_iterator
+ (&triangulation,
+ 1,
+ 4,
+ &dof_handler));
+
+ std::vector<Vector<double> > m (points.size());
+ fe_function.vector_value_list (points, m);
+
+ for (unsigned int i=0; i<m.size(); ++i)
+ {
+ Assert (std::fabs(m[i](0) - points[i].square())
+ <
+ 1e-10 * std::fabs(m[i](0) + points[i].square()),
+ ExcInternalError());
+
+ Assert (std::fabs(m[i](1))
+ <
+ 1e-10,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+int main ()
+{
+ std::ofstream logfile("fe_field_function_05_vector/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ test<2>();
+
+ return 0;
+}
+
--- /dev/null
+
+DEAL::OK
projects / dealii-svn.git / commitdiff