anna_1.exe : anna_1.g.$(OBJEXT) $(libraries)
anna_2.exe : anna_2.g.$(OBJEXT) $(libraries)
anna_3.exe : anna_3.g.$(OBJEXT) $(libraries)
+anna_4.exe : anna_4.g.$(OBJEXT) $(libraries)
anna_5.exe : anna_5.g.$(OBJEXT) $(libraries)
gerold_1.exe : gerold_1.g.$(OBJEXT) $(libraries)
roy_1.exe : roy_1.g.$(OBJEXT) $(libraries)
denis_1.exe : denis_1.g.$(OBJEXT) $(libraries)
+unit_support_points.exe : unit_support_points.g.$(OBJEXT) $(libraries)
-tests = anna_1 anna_2 anna_3 anna_5 \
+tests = anna_1 anna_2 anna_3 anna_4 anna_5 \
gerold_1 \
roy_1 \
- denis_1
+ denis_1 \
+ unit_support_points
############################################################
--- /dev/null
+//---------------------------- anna_4.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002 by the deal.II authors and Anna Schneebeli
+//
+// 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.
+//
+//---------------------------- anna_4.cc ---------------------------
+
+
+// check for an abort in VectorTools::interpolate_boundary_values.
+//
+// this program is a modified version of one by Anna Schneebeli,
+// University of Basel
+
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+
+#include <grid/tria.h>
+#include <dofs/dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <numerics/vectors.h>
+
+ // We need a FESystem
+#include <fe/fe_system.h>
+
+ // we need DG-elements
+ // and Q1-elements
+#include <fe/fe_q.h>
+#include <fe/fe_dgp.h>
+
+#include <fstream>
+#include <iostream>
+#include <vector>
+
+
+template <int dim>
+class VectorBoundaryValues : public Function<dim>
+{
+ public:
+ VectorBoundaryValues ();
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &values) const;
+};
+
+template <int dim>
+VectorBoundaryValues<dim>::VectorBoundaryValues () :
+ Function<dim> (2)
+{};
+
+template <int dim>
+inline
+void VectorBoundaryValues<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+{
+ Assert (values.size() == 2,
+ ExcDimensionMismatch (values.size(), 2));
+
+ for (unsigned int i=0; i<2; ++i)
+ values(i) = p(i)*p(i);
+};
+
+
+
+template <int dim>
+class FindBug
+{
+ public:
+ FindBug ();
+ void run ();
+ private:
+ void make_grid_and_dofs ();
+ void dirichlet_conditions ();
+
+ Triangulation<dim> triangulation;
+ FESystem<dim> fe;
+ DoFHandler<dim> dof_handler;
+ Vector<double> solution;
+};
+
+
+ // Construct FESystem with
+ // first component: Q1-Element,
+ // second component: lowest order DG_Element
+template <int dim>
+FindBug<dim>::FindBug () :
+ fe (FE_Q<dim>(1), 1,
+ FE_DGP<dim>(0), 1),
+ dof_handler (triangulation)
+{};
+
+
+template <int dim>
+void FindBug<dim>::make_grid_and_dofs ()
+{
+
+ GridGenerator::hyper_cube (triangulation);
+ triangulation.refine_global (1);
+
+ deallog << "Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl;
+
+ deallog << "Total number of cells: "
+ << triangulation.n_cells()
+ << std::endl;
+
+
+ dof_handler.distribute_dofs (fe);
+
+
+ deallog << "Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ solution.reinit(dof_handler.n_dofs());
+
+};
+
+
+template <int dim>
+void FindBug<dim>::dirichlet_conditions ()
+{
+ // we want to set the Boundary DoFs
+ // of the selected component to a
+ // given value, say zero. To do
+ // so, we want to use VectorTools::
+ // interpolate_boundary_values
+ //
+ // This works fine if all the
+ // components have support on the
+ // faces. (This, of course, has to
+ // be requested when fixing the
+ // boundary DoFs.) However,
+ // getting the values for the
+ // boundary DoFs of a valid
+ // component by the function
+ // VectorTools::
+ // interpolate_boundary_values and
+ // an appropriate component_mask
+ // (unselecting the DG-component)
+ // does not work yet. Or, better:
+ // it did not work and this program
+ // checks that it has been
+ // correctly implemented by now
+
+
+ std::map<unsigned int,double> dirichlet_dofs;
+
+ // we declare a vector of bools,
+ // which tells the
+ // VectorTools::interpolate_boundary_values
+ // on which components of the
+ // FESystem we want to impose
+ // Dirichlet BC.
+ std::vector<bool> component_mask(2);
+ // Dirichlet-BC for the
+ // Q1-Component
+ component_mask[0] = true;
+ // no Dirichlet-BC for the
+ // DG-component
+ component_mask[1] = false;
+
+ // This is just for the final
+ // output-test
+ for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+ dirichlet_dofs[i] = 1.;
+
+
+ // Here comes the crucial call....
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim> (2),
+ dirichlet_dofs,
+ component_mask);
+
+
+ std::vector<bool> fixed_dofs (dof_handler.n_dofs());
+ std::set<unsigned char> boundary_indicators;
+ boundary_indicators.insert (0);
+
+ // get a list of those boundary DoFs which
+ // we want to be fixed:
+ DoFTools::extract_boundary_dofs (dof_handler,
+ component_mask,
+ fixed_dofs,
+ boundary_indicators);
+
+ // (Primitive) Check if the DoFs
+ // where adjusted correctly (note
+ // that we have preset all values
+ // to 1, and interpolate_b_v should
+ // have overwritten those for
+ // component 0 by 0)
+ for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+ {
+ if (fixed_dofs[i] == true)
+ {
+ Assert (dirichlet_dofs[i] == 0, ExcInternalError());
+ }
+ else
+ {
+ Assert (dirichlet_dofs[i] == 1, ExcInternalError());
+ };
+ };
+
+ // check 1 has obviously succeeded,
+ // so also check a more complicated
+ // boundary value function
+ dirichlet_dofs.clear ();
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ VectorBoundaryValues<dim> (),
+ dirichlet_dofs,
+ component_mask);
+ for (unsigned int i=0; i<dof_handler.n_dofs(); ++i)
+ if (fixed_dofs[i] == true)
+ deallog << i << " " << dirichlet_dofs[i] << std::endl;
+};
+
+
+
+template <int dim>
+void FindBug<dim>::run ()
+{
+ make_grid_and_dofs ();
+ dirichlet_conditions ();
+};
+
+
+
+int main ()
+{
+ std::ofstream logfile("anna_4.output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ FindBug<2>().run ();
+ FindBug<3>().run ();
+
+ return 0;
+};
+
--- /dev/null
+
+DEAL::Number of active cells: 4
+DEAL::Total number of cells: 5
+DEAL::Number of degrees of freedom: 13
+DEAL::0 0.00000
+DEAL::1 0.250000
+DEAL::3 0.00000
+DEAL::5 1.00000
+DEAL::6 1.00000
+DEAL::8 1.00000
+DEAL::9 0.250000
+DEAL::11 0.00000
+DEAL::Number of active cells: 8
+DEAL::Total number of cells: 9
+DEAL::Number of degrees of freedom: 35
+DEAL::0 0.00000
+DEAL::1 0.250000
+DEAL::2 0.250000
+DEAL::3 0.00000
+DEAL::4 0.00000
+DEAL::5 0.250000
+DEAL::7 0.00000
+DEAL::9 1.00000
+DEAL::10 1.00000
+DEAL::11 1.00000
+DEAL::12 1.00000
+DEAL::14 1.00000
+DEAL::15 0.250000
+DEAL::16 1.00000
+DEAL::17 0.250000
+DEAL::19 0.00000
+DEAL::20 0.00000
+DEAL::22 0.00000
+DEAL::23 0.250000
+DEAL::24 0.250000
+DEAL::25 0.00000
+DEAL::27 1.00000
+DEAL::28 1.00000
+DEAL::30 1.00000
+DEAL::31 0.250000
+DEAL::33 0.00000
--- /dev/null
+//---------------------------- unit_support_points.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2002 by the deal.II authors and Anna Schneebeli
+//
+// 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.
+//
+//---------------------------- unit_support_points.cc ---------------------------
+
+
+// check equivalence of get_unit_support_points() and
+// unit_support_point(), as well as the case where an element does not
+// have support points
+
+#include <base/logstream.h>
+#include <fe/fe_system.h>
+#include <fe/fe_q.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_dgp.h>
+#include <fe/fe_nedelec.h>
+#include <fstream>
+
+
+template <int dim>
+void check_cell1 (const FiniteElement<dim> &fe)
+{
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ Assert (fe.get_unit_support_points()[i] ==
+ fe.unit_support_point(i),
+ ExcInternalError());
+ deallog << "dim=" << dim << ", cell=ok" << std::endl;
+};
+
+
+template <int dim>
+void check_face1 (const FiniteElement<dim> &fe)
+{
+ for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+ Assert (fe.get_unit_face_support_points()[i] ==
+ fe.unit_face_support_point(i),
+ ExcInternalError());
+ deallog << "dim=" << dim << ", face=ok" << std::endl;
+};
+
+
+void check_face1 (const FiniteElement<1> &)
+{};
+
+
+
+template <int dim>
+void check1 (const FiniteElement<dim> &fe)
+{
+ check_cell1 (fe);
+ check_face1 (fe);
+};
+
+
+template <int dim>
+void check_cell2 (const FiniteElement<dim> &fe,
+ const unsigned int comp)
+{
+ for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+ if (fe.system_to_component_index(i).first == comp)
+ deallog << i << " " << fe.unit_support_point(i)
+ << std::endl;
+ deallog << "dim=" << dim << ", cell=ok" << std::endl;
+};
+
+
+template <int dim>
+void check_face2 (const FiniteElement<dim> &fe,
+ const unsigned int comp)
+{
+ for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+ if (fe.face_system_to_component_index(i).first == comp)
+ deallog << i << " " << fe.unit_face_support_point(i)
+ << std::endl;
+ deallog << "dim=" << dim << ", face=ok" << std::endl;
+};
+
+
+void check_face2 (const FiniteElement<1> &,
+ const unsigned int)
+{};
+
+
+
+template <int dim>
+void check2 (const FiniteElement<dim> &fe,
+ const unsigned int comp)
+{
+ check_cell2 (fe, comp);
+ check_face2 (fe, comp);
+};
+
+
+
+template <int dim>
+void check ()
+{
+ check1 (FE_Q<dim>(2));
+ check1 (FE_DGQ<dim>(2));
+ check1 (FESystem<dim> (FE_Q<dim> (2), 2,
+ FE_DGQ<dim> (2),1));
+
+ check1 (FESystem<dim> (FE_Q<dim> (2), 2,
+ FE_DGQ<dim> (2),1));
+
+ check2 (FESystem<dim> (FE_Q<dim> (2), 1,
+ FE_DGQ<dim> (2), 1),
+ 1);
+ check2 (FESystem<dim> (FE_Q<dim> (2), 1,
+ FE_DGP<dim> (2), 1),
+ 0);
+};
+
+
+
+int main ()
+{
+ std::ofstream logfile("unit_support_points.output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ check<1> ();
+ check<2> ();
+ check<3> ();
+ return 0;
+};
--- /dev/null
+
+DEAL::dim=1, cell=ok
+DEAL::dim=1, cell=ok
+DEAL::dim=1, cell=ok
+DEAL::dim=1, cell=ok
+DEAL::3 0.00000
+DEAL::4 0.500000
+DEAL::5 1.00000
+DEAL::dim=1, cell=ok
+DEAL::0 0.00000
+DEAL::1 1.00000
+DEAL::2 0.500000
+DEAL::dim=1, cell=ok
+DEAL::dim=2, cell=ok
+DEAL::dim=2, face=ok
+DEAL::dim=2, cell=ok
+DEAL::dim=2, face=ok
+DEAL::dim=2, cell=ok
+DEAL::dim=2, face=ok
+DEAL::dim=2, cell=ok
+DEAL::dim=2, face=ok
+DEAL::9 0.00000 0.00000
+DEAL::10 0.500000 0.00000
+DEAL::11 1.00000 0.00000
+DEAL::12 0.00000 0.500000
+DEAL::13 0.500000 0.500000
+DEAL::14 1.00000 0.500000
+DEAL::15 0.00000 1.00000
+DEAL::16 0.500000 1.00000
+DEAL::17 1.00000 1.00000
+DEAL::dim=2, cell=ok
+DEAL::dim=2, face=ok
+DEAL::0 0.00000 0.00000
+DEAL::1 1.00000 0.00000
+DEAL::2 1.00000 1.00000
+DEAL::3 0.00000 1.00000
+DEAL::4 0.500000 0.00000
+DEAL::5 1.00000 0.500000
+DEAL::6 0.500000 1.00000
+DEAL::7 0.00000 0.500000
+DEAL::8 0.500000 0.500000
+DEAL::dim=2, cell=ok
+DEAL::0 0.00000
+DEAL::1 1.00000
+DEAL::2 0.500000
+DEAL::dim=2, face=ok
+DEAL::dim=3, cell=ok
+DEAL::dim=3, face=ok
+DEAL::dim=3, cell=ok
+DEAL::dim=3, face=ok
+DEAL::dim=3, cell=ok
+DEAL::dim=3, face=ok
+DEAL::dim=3, cell=ok
+DEAL::dim=3, face=ok
+DEAL::27 0.00000 0.00000 0.00000
+DEAL::28 0.500000 0.00000 0.00000
+DEAL::29 1.00000 0.00000 0.00000
+DEAL::30 0.00000 0.500000 0.00000
+DEAL::31 0.500000 0.500000 0.00000
+DEAL::32 1.00000 0.500000 0.00000
+DEAL::33 0.00000 1.00000 0.00000
+DEAL::34 0.500000 1.00000 0.00000
+DEAL::35 1.00000 1.00000 0.00000
+DEAL::36 0.00000 0.00000 0.500000
+DEAL::37 0.500000 0.00000 0.500000
+DEAL::38 1.00000 0.00000 0.500000
+DEAL::39 0.00000 0.500000 0.500000
+DEAL::40 0.500000 0.500000 0.500000
+DEAL::41 1.00000 0.500000 0.500000
+DEAL::42 0.00000 1.00000 0.500000
+DEAL::43 0.500000 1.00000 0.500000
+DEAL::44 1.00000 1.00000 0.500000
+DEAL::45 0.00000 0.00000 1.00000
+DEAL::46 0.500000 0.00000 1.00000
+DEAL::47 1.00000 0.00000 1.00000
+DEAL::48 0.00000 0.500000 1.00000
+DEAL::49 0.500000 0.500000 1.00000
+DEAL::50 1.00000 0.500000 1.00000
+DEAL::51 0.00000 1.00000 1.00000
+DEAL::52 0.500000 1.00000 1.00000
+DEAL::53 1.00000 1.00000 1.00000
+DEAL::dim=3, cell=ok
+DEAL::dim=3, face=ok
+DEAL::0 0.00000 0.00000 0.00000
+DEAL::1 1.00000 0.00000 0.00000
+DEAL::2 1.00000 0.00000 1.00000
+DEAL::3 0.00000 0.00000 1.00000
+DEAL::4 0.00000 1.00000 0.00000
+DEAL::5 1.00000 1.00000 0.00000
+DEAL::6 1.00000 1.00000 1.00000
+DEAL::7 0.00000 1.00000 1.00000
+DEAL::8 0.500000 0.00000 0.00000
+DEAL::9 1.00000 0.00000 0.500000
+DEAL::10 0.500000 0.00000 1.00000
+DEAL::11 0.00000 0.00000 0.500000
+DEAL::12 0.500000 1.00000 0.00000
+DEAL::13 1.00000 1.00000 0.500000
+DEAL::14 0.500000 1.00000 1.00000
+DEAL::15 0.00000 1.00000 0.500000
+DEAL::16 0.00000 0.500000 0.00000
+DEAL::17 1.00000 0.500000 0.00000
+DEAL::18 1.00000 0.500000 1.00000
+DEAL::19 0.00000 0.500000 1.00000
+DEAL::20 0.500000 0.00000 0.500000
+DEAL::21 0.500000 1.00000 0.500000
+DEAL::22 0.500000 0.500000 0.00000
+DEAL::23 1.00000 0.500000 0.500000
+DEAL::24 0.500000 0.500000 1.00000
+DEAL::25 0.00000 0.500000 0.500000
+DEAL::26 0.500000 0.500000 0.500000
+DEAL::dim=3, cell=ok
+DEAL::0 0.00000 0.00000
+DEAL::1 1.00000 0.00000
+DEAL::2 1.00000 1.00000
+DEAL::3 0.00000 1.00000
+DEAL::4 0.500000 0.00000
+DEAL::5 1.00000 0.500000
+DEAL::6 0.500000 1.00000
+DEAL::7 0.00000 0.500000
+DEAL::8 0.500000 0.500000
+DEAL::dim=3, face=ok