mesh_3d_10.exe : mesh_3d_10.g.$(OBJEXT) $(libraries)
mesh_3d_11.exe : mesh_3d_11.g.$(OBJEXT) $(libraries)
mesh_3d_12.exe : mesh_3d_12.g.$(OBJEXT) $(libraries)
+mesh_3d_13.exe : mesh_3d_13.g.$(OBJEXT) $(libraries)
+mesh_3d_14.exe : mesh_3d_14.g.$(OBJEXT) $(libraries)
q_points.exe : q_points.g.$(OBJEXT) $(libraries)
hyper_ball_3d cylinder coarsening_3d \
mesh_3d_1 mesh_3d_2 mesh_3d_3 mesh_3d_4 mesh_3d_5 \
mesh_3d_6 mesh_3d_7 mesh_3d_8 mesh_3d_9 mesh_3d_10 \
- mesh_3d_11 mesh_3d_12 \
+ mesh_3d_11 mesh_3d_12 mesh_3d_13 mesh_3d_14 \
q_points
############################################################
// the Kelly estimator used to fall over when presented with 3d meshes
-// with mis-oriented faces
+// with mis-oriented faces. the actual assertion where we fail is
+// again tested isolated in mesh_3d_13 (we limit the number of
+// refinements here to reduce run time; it is not limited there to the
+// same degree)
#include "mesh_3d.h"
#include <base/logstream.h>
#include <base/quadrature_lib.h>
+#include <base/function_lib.h>
#include <lac/vector.h>
#include <grid/tria.h>
#include <grid/tria_accessor.h>
#include <dofs/dof_handler.h>
#include <fe/fe_q.h>
#include <numerics/error_estimator.h>
+#include <numerics/vectors.h>
#include <fstream>
+
+
void check_this (Triangulation<3> &tria)
{
FE_Q<3> fe(1);
Vector<double> u(dof_handler.n_dofs());
Vector<float> e(tria.n_active_cells());
+
+ VectorTools::interpolate (dof_handler,
+ Functions::SquareFunction<3>(),
+ u);
KellyErrorEstimator<3>::estimate (dof_handler,
QGauss2<2>(),
FunctionMap<3>::type(),
u,
e);
+
+ deallog << " " << e.l1_norm() << std::endl;
+ deallog << " " << e.l2_norm() << std::endl;
+ deallog << " " << e.linfty_norm() << std::endl;
}
deallog << "Initial check" << std::endl;
check_this (tria);
- for (unsigned int r=0; r<3; ++r)
+ for (unsigned int r=0; r<2; ++r)
{
tria.refine_global (1);
deallog << "Check " << r << std::endl;
check_this (tria);
}
-
- coarsen_global (tria);
- deallog << "Check " << 1 << std::endl;
- check_this (tria);
-
- tria.refine_global (1);
- deallog << "Check " << 2 << std::endl;
- check_this (tria);
}
--- /dev/null
+//---------------------------- mesh_3d_13.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2003 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.
+//
+//---------------------------- mesh_3d_13.cc ---------------------------
+
+
+// this tests the assertion that is actually triggered in mesh_3d_12,
+// isolated from the rest of the code in which it sits in that test
+//
+// actually, the computation of the orientation flag was wrong as we
+// were also considering the orientation flag of the present cell,
+// which we shouldn't have
+
+#include "mesh_3d.h"
+
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_reordering.h>
+#include <grid/grid_generator.h>
+
+#include <fstream>
+
+
+void check_this (Triangulation<3> &tria)
+{
+ Triangulation<3>::active_cell_iterator cell = tria.begin_active();
+ for (; cell!=tria.end(); ++cell)
+ for (unsigned int face_no=0; face_no<GeometryInfo<3>::faces_per_cell;
+ ++face_no)
+ if (!cell->at_boundary(face_no)
+ &&
+ cell->neighbor(face_no)->has_children())
+ for (unsigned int subface_no=0;
+ subface_no<GeometryInfo<3>::subfaces_per_face;
+ ++subface_no)
+ {
+ // get an iterator
+ // pointing to the cell
+ // behind the present
+ // subface
+ const Triangulation<3>::cell_iterator
+ neighbor = cell->neighbor(face_no);
+ const unsigned int neighbor_neighbor
+ = cell->neighbor_of_neighbor (face_no);
+ const bool orientation_flag
+ = (neighbor->face_orientation(neighbor_neighbor) == true);
+ static const unsigned int subface_translation[4]
+ = { 0, 3, 2, 1 };
+ const unsigned int neighbor_child_index
+ = (GeometryInfo<3>::
+ child_cell_on_face(neighbor_neighbor,
+ (orientation_flag ?
+ subface_no :
+ subface_translation[subface_no])));
+ const Triangulation<3>::active_cell_iterator neighbor_child
+ = neighbor->child(neighbor_child_index);
+
+ Assert (neighbor_child->face(neighbor_neighbor) ==
+ cell->face(face_no)->child(subface_no),
+ ExcInternalError());
+ Assert (!neighbor->child(neighbor_child_index)->has_children(),
+ ExcInternalError());
+ }
+}
+
+
+
+void check (Triangulation<3> &tria)
+{
+ (++tria.begin_active())->set_refine_flag ();
+ tria.execute_coarsening_and_refinement ();
+
+ deallog << "Initial check" << std::endl;
+ check_this (tria);
+
+ for (unsigned int r=0; r<3; ++r)
+ {
+ tria.refine_global (1);
+ deallog << "Check " << r << std::endl;
+ check_this (tria);
+ }
+
+ coarsen_global (tria);
+ deallog << "Check " << 1 << std::endl;
+ check_this (tria);
+
+ tria.refine_global (1);
+ deallog << "Check " << 2 << std::endl;
+ check_this (tria);
+}
+
+
+int main ()
+{
+ std::ofstream logfile("mesh_3d_13.output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ {
+ Triangulation<3> coarse_grid;
+ create_two_cubes (coarse_grid);
+ check (coarse_grid);
+ }
+
+ {
+ Triangulation<3> coarse_grid;
+ create_L_shape (coarse_grid);
+ check (coarse_grid);
+ }
+
+ {
+ Triangulation<3> coarse_grid;
+ GridGenerator::hyper_ball (coarse_grid);
+ check (coarse_grid);
+ }
+
+}
+
+
+
--- /dev/null
+//---------------------------- mesh_3d_14.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2003 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.
+//
+//---------------------------- mesh_3d_14.cc ---------------------------
+
+
+// similar to mesh_3d_13, but check that quadrature points are
+// correct. it is thus also similar to mesh_3d_7, but checks for faces
+// and subfaces
+
+#include "mesh_3d.h"
+
+#include <base/logstream.h>
+#include <base/quadrature_lib.h>
+#include <lac/vector.h>
+#include <grid/tria.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/grid_reordering.h>
+#include <grid/grid_generator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <fe/fe_values.h>
+
+#include <fstream>
+
+
+void check_this (Triangulation<3> &tria)
+{
+ QTrapez<2> quadrature;
+ FE_Q<3> fe(1);
+ FEFaceValues<3> fe_face_values1 (fe, quadrature,
+ update_q_points | update_JxW_values |
+ update_normal_vectors);
+ FESubfaceValues<3> fe_face_values2 (fe, quadrature,
+ update_q_points | update_JxW_values |
+ update_normal_vectors);
+
+ DoFHandler<3> dof_handler (tria);
+ dof_handler.distribute_dofs (fe);
+
+ DoFHandler<3>::active_cell_iterator cell = dof_handler.begin_active();
+ for (; cell!=dof_handler.end(); ++cell)
+ for (unsigned int face_no=0; face_no<GeometryInfo<3>::faces_per_cell;
+ ++face_no)
+ if (!cell->at_boundary(face_no)
+ &&
+ cell->neighbor(face_no)->has_children())
+ for (unsigned int subface_no=0;
+ subface_no<GeometryInfo<3>::subfaces_per_face;
+ ++subface_no)
+ {
+ // get an iterator
+ // pointing to the cell
+ // behind the present
+ // subface
+ const DoFHandler<3>::cell_iterator
+ neighbor = cell->neighbor(face_no);
+ const unsigned int neighbor_neighbor
+ = cell->neighbor_of_neighbor (face_no);
+
+ // see whether face and
+ // the neighbor's
+ // counterface share the
+ // same indexing of
+ // children. if not so,
+ // translate child
+ // indices
+ const bool face_orientations_match
+ = (neighbor->face_orientation(neighbor_neighbor) ==
+ cell->face_orientation(face_no));
+ static const unsigned int subface_translation[4]
+ = { 0, 3, 2, 1 };
+ const unsigned int neighbor_child_index
+ = (GeometryInfo<3>::
+ child_cell_on_face(neighbor_neighbor,
+ (face_orientations_match ?
+ subface_no :
+ subface_translation[subface_no])));
+ const DoFHandler<3>::active_cell_iterator neighbor_child
+ = neighbor->child(neighbor_child_index);
+
+ fe_face_values1.reinit (neighbor_child, neighbor_neighbor);
+ fe_face_values2.reinit (cell, face_no, subface_no);
+
+ for (unsigned int q=0; q<quadrature.n_quadrature_points; ++q)
+ {
+ Assert ((fe_face_values1.quadrature_point(q)-
+ fe_face_values2.quadrature_point(q)).square()
+ < 1e-20,
+ ExcInternalError());
+
+ Assert (std::fabs(fe_face_values1.JxW(q)-
+ fe_face_values2.JxW(q)) < 1e-15,
+ ExcInternalError());
+ Assert ((fe_face_values1.normal_vector(q) +
+ fe_face_values2.normal_vector(q)).square()
+ < 1e-20,
+ ExcInternalError());
+ }
+ }
+}
+
+
+
+void check (Triangulation<3> &tria)
+{
+ (++tria.begin_active())->set_refine_flag ();
+ tria.execute_coarsening_and_refinement ();
+
+ deallog << "Initial check" << std::endl;
+ check_this (tria);
+
+ for (unsigned int r=0; r<3; ++r)
+ {
+ tria.refine_global (1);
+ deallog << "Check " << r << std::endl;
+ check_this (tria);
+ }
+
+ coarsen_global (tria);
+ deallog << "Check " << 1 << std::endl;
+ check_this (tria);
+
+ tria.refine_global (1);
+ deallog << "Check " << 2 << std::endl;
+ check_this (tria);
+}
+
+
+int main ()
+{
+ std::ofstream logfile("mesh_3d_14.output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ {
+ Triangulation<3> coarse_grid;
+ create_two_cubes (coarse_grid);
+ check (coarse_grid);
+ }
+
+ {
+ Triangulation<3> coarse_grid;
+ create_L_shape (coarse_grid);
+ check (coarse_grid);
+ }
+
+ {
+ Triangulation<3> coarse_grid;
+ GridGenerator::hyper_ball (coarse_grid);
+ check (coarse_grid);
+ }
+
+}
+
+
+
--- /dev/null
+
+DEAL::Initial check
+DEAL:: 1.93150
+DEAL:: 0.678289
+DEAL:: 0.402964
+DEAL::Check 0
+DEAL:: 3.77915
+DEAL:: 0.572837
+DEAL:: 0.179270
+DEAL::Check 1
+DEAL:: 5.86911
+DEAL:: 0.329660
+DEAL:: 0.0411273
+DEAL::Initial check
+DEAL:: 2.52567
+DEAL:: 0.838836
+DEAL:: 0.402964
+DEAL::Check 0
+DEAL:: 5.28305
+DEAL:: 0.758350
+DEAL:: 0.179270
+DEAL::Check 1
+DEAL:: 8.26356
+DEAL:: 0.440081
+DEAL:: 0.0411273
+DEAL::Initial check
+DEAL:: 6.63199
+DEAL:: 2.18605
+DEAL:: 0.971172
+DEAL::Check 0
+DEAL:: 7.24766
+DEAL:: 0.921457
+DEAL:: 0.184524
+DEAL::Check 1
+DEAL:: 9.04235
+DEAL:: 0.421062
+DEAL:: 0.0400020
--- /dev/null
+
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
--- /dev/null
+
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
--- /dev/null
+
+DEAL::Initial check
+DEAL:: 1.93150
+DEAL:: 0.678289
+DEAL:: 0.402964
+DEAL::Check 0
+DEAL:: 3.77915
+DEAL:: 0.572837
+DEAL:: 0.179270
+DEAL::Check 1
+DEAL:: 5.86911
+DEAL:: 0.329660
+DEAL:: 0.0411273
+DEAL::Initial check
+DEAL:: 2.52567
+DEAL:: 0.838836
+DEAL:: 0.402964
+DEAL::Check 0
+DEAL:: 5.28305
+DEAL:: 0.758350
+DEAL:: 0.179270
+DEAL::Check 1
+DEAL:: 8.26356
+DEAL:: 0.440081
+DEAL:: 0.0411273
+DEAL::Initial check
+DEAL:: 6.63199
+DEAL:: 2.18605
+DEAL:: 0.971172
+DEAL::Check 0
+DEAL:: 7.24766
+DEAL:: 0.921457
+DEAL:: 0.184524
+DEAL::Check 1
+DEAL:: 9.04235
+DEAL:: 0.421062
+DEAL:: 0.0400020
--- /dev/null
+
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
--- /dev/null
+
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
--- /dev/null
+
+DEAL::Initial check
+DEAL:: 1.93150
+DEAL:: 0.678289
+DEAL:: 0.402964
+DEAL::Check 0
+DEAL:: 3.77915
+DEAL:: 0.572837
+DEAL:: 0.179270
+DEAL::Check 1
+DEAL:: 5.86911
+DEAL:: 0.329660
+DEAL:: 0.0411273
+DEAL::Initial check
+DEAL:: 2.52567
+DEAL:: 0.838836
+DEAL:: 0.402964
+DEAL::Check 0
+DEAL:: 5.28305
+DEAL:: 0.758350
+DEAL:: 0.179270
+DEAL::Check 1
+DEAL:: 8.26356
+DEAL:: 0.440081
+DEAL:: 0.0411273
+DEAL::Initial check
+DEAL:: 6.63199
+DEAL:: 2.18605
+DEAL:: 0.971172
+DEAL::Check 0
+DEAL:: 7.24766
+DEAL:: 0.921457
+DEAL:: 0.184524
+DEAL::Check 1
+DEAL:: 9.04235
+DEAL:: 0.421062
+DEAL:: 0.0400020
--- /dev/null
+
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
--- /dev/null
+
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Initial check
+DEAL::Check 0
+DEAL::Check 1
+DEAL::Check 2
+DEAL::Check 1
+DEAL::Check 2