From f63dac72543d4fb8580a38f929732d2355dc8a0e Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 13 Jan 2015 14:16:41 -0600 Subject: [PATCH] Remove old benchmarks that haven't been looked at in years. They almost certainly also don't build any more. --- tests/benchmarks/old/Makefile | 83 - tests/benchmarks/old/benchmark.cc | 72 - tests/benchmarks/old/dof_handler_size.cc | 120 - tests/benchmarks/old/dof_handler_timing.cc | 215 - tests/benchmarks/old/dof_handler_timing_01.cc | 116 - .../benchmarks/old/patches/library/2007-05-30 | 91 - tests/benchmarks/old/performance.cc | 141 - tests/benchmarks/old/performance.pl | 27 - tests/benchmarks/old/spec2006-447.dealII.cc | 4279 ----------------- 9 files changed, 5144 deletions(-) delete mode 100644 tests/benchmarks/old/Makefile delete mode 100644 tests/benchmarks/old/benchmark.cc delete mode 100644 tests/benchmarks/old/dof_handler_size.cc delete mode 100644 tests/benchmarks/old/dof_handler_timing.cc delete mode 100644 tests/benchmarks/old/dof_handler_timing_01.cc delete mode 100644 tests/benchmarks/old/patches/library/2007-05-30 delete mode 100644 tests/benchmarks/old/performance.cc delete mode 100644 tests/benchmarks/old/performance.pl delete mode 100644 tests/benchmarks/old/spec2006-447.dealII.cc diff --git a/tests/benchmarks/old/Makefile b/tests/benchmarks/old/Makefile deleted file mode 100644 index 125462fc7e..0000000000 --- a/tests/benchmarks/old/Makefile +++ /dev/null @@ -1,83 +0,0 @@ -############################################################ -# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors -############################################################ - -############################################################ -# Include general settings for including DEAL libraries -############################################################ - -include ../Makefile.paths -include $D/common/Make.global_options -debug-mode = on - -libraries = $(lib-deal2-2d.o) $(lib-deal2-3d.o) $(lib-lac.o) $(lib-base.o) -libraries.g = $(lib-deal2-2d.g) $(lib-deal2-3d.g) $(lib-lac.g) $(lib-base.g) - -default: run-tests - -############################################################ - -tests_x = dof_handler_* - -# from above list of regular expressions, generate the real set of -# tests -expand = $(shell echo $(addsuffix .cc,$(1)) \ - | $(PERL) -pi -e 's/\.cc//g;') -tests = $(call expand,$(tests_x)) - - -############################################################ -# First how to create executables, including all necessary -# flags: -############################################################ - -flags = $(CXXFLAGS.g) - -ifeq ($(findstring gcc,$(GXX_VERSION)),gcc) -flags += -Wno-missing-noreturn -endif - -%/obj.g.$(OBJEXT) : %.cc - @echo =====debug========= $< - @echo Compiling > $*/status - @$(CXX) $(flags) -c $< -o $@ - -%/obj.$(OBJEXT) : %.cc - @echo =====optimized===== $< - @echo Compiling > $*/status - @$(CXX) $(CXXFLAGS.o) -c $< -o $@ - -%/exe : %/obj.$(OBJEXT) $(libraries) - @echo =====linking======= $@ - @echo Linking > $*/status - @$(CXX) $(LDFLAGS) -o $@ $^ $(LIBS) - -%/g.exe : %/obj.g.$(OBJEXT) $(libraries.g) - @echo =====linking======= $@ - @echo Linking > $*/status - @$(CXX) $(LDFLAGS) -o $@ $^ $(LIBS) - -%/output : %/exe - @echo =====running======= $< - @$< - -############################################################ -# After all these general rules, here is the target to be -# executed by make: for each entry in the list $(tests) -# perform a check. -############################################################ -run-tests: $(tests:%=%/output) - -build: $(tests:%=%/exe) - -refs: $(tests:%=%/ref) - - -Makefile.depend: $(shell echo *.cc) - @echo =====Dependencies== $@ - @$D/common/scripts/make_dependencies $(INCLUDE) $^ > $@ \ - || (rm -f $@ ; false) - @$(PERL) -pi -e 's#\.((g\.)?$(OBJEXT):)#/obj.\1#g;' $@ - - --include Makefile.depend diff --git a/tests/benchmarks/old/benchmark.cc b/tests/benchmarks/old/benchmark.cc deleted file mode 100644 index aa403c4e20..0000000000 --- a/tests/benchmarks/old/benchmark.cc +++ /dev/null @@ -1,72 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 1998 - 2013 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - - -#include "../tests.h" -#include -#include -#include "quickmatrix.h" -#include - -#define ITER 100 - - -int main() -{ - Vector u; - Vector v; - - clock_t start; - clock_t diff; - - deallog << "Iterations: " << ITER << std::endl; - - for (unsigned int nx=32; nx<8192 ; nx*=2) - { - const unsigned int dim=(nx-1)*(nx-1); - - deallog << "size = " << nx << " dim = " << dim << std::endl; - - start = clock(); - for (unsigned int i=0; i A(nx,nx); - - start = clock(); - for (unsigned int i=0; i -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace dealii; - - -// Fill dof handlers for different elements and see how large they get. -template -void check_dofs(DOF &dof) -{ - FE_Q q1(1); - FE_Q q3(3); - FESystem sys1(q3, 1); - FESystem sys2(q3, 10); - - dof.distribute_dofs(q1); - deallog << "Dofs Q1 " << std::setw(12) << dof.n_dofs() - << " memory " << std::setw(12) << dof.memory_consumption() - << " quotient " << (1./dof.n_dofs()*dof.memory_consumption()) - << std::endl; - - dof.distribute_dofs(q3); - deallog << "Dofs Q3 " << std::setw(12) << dof.n_dofs() - << " memory " << std::setw(12) << dof.memory_consumption() - << " quotient " << (1./dof.n_dofs()*dof.memory_consumption()) - << std::endl; - - dof.distribute_dofs(sys1); - deallog << "Dofs Sys1 " << std::setw(12) << dof.n_dofs() - << " memory " << std::setw(12) << dof.memory_consumption() - << " quotient " << (1./dof.n_dofs()*dof.memory_consumption()) - << std::endl; - - dof.distribute_dofs(sys2); - deallog << "Dofs Sys2 " << std::setw(12) << dof.n_dofs() - << " memory " << std::setw(12) << dof.memory_consumption() - << " quotient " << (1./dof.n_dofs()*dof.memory_consumption()) - << std::endl; - - dof.clear(); -} - - -template -void check (bool local) -{ - deallog << "Dimension " << dim << std::endl; - Triangulation tr(Triangulation::maximum_smoothing); - GridGenerator::hyper_cube(tr); - - if (local) - for (unsigned int i=0; i<99/dim; ++i) - { - tr.last_active()->set_refine_flag(); - tr.execute_coarsening_and_refinement(); - } - else - tr.refine_global(18/dim); - - deallog << "Levels " << tr.n_levels() << "Cells "<< tr.n_cells() - << std::endl - << " active " << std::setw(12) << tr.n_active_cells() - << " memory " << std::setw(12) << tr.memory_consumption() - << " quotient " << (1./tr.n_cells()*tr.memory_consumption()) - << std::endl; - - FE_Q q1(1); - FE_Q q3(3); - FESystem sys1(q3, 1); - FESystem sys2(q3, 10); - deallog.push("DoF"); - DoFHandler dof(tr); - check_dofs(dof); - deallog.pop(); - deallog.push("MGDoF"); - DoFHandler mgdof(tr); - check_dofs(mgdof); - deallog.pop(); -} - -int main() -{ - std::ofstream out("dof_handler_size/output"); - deallog.attach(out); - deallog.push("local"); - check<2>(true); - check<3>(true); - deallog.pop(); - deallog.push("global"); - check<2>(false); - check<3>(false); - deallog.pop(); -} - diff --git a/tests/benchmarks/old/dof_handler_timing.cc b/tests/benchmarks/old/dof_handler_timing.cc deleted file mode 100644 index 87e3883ce1..0000000000 --- a/tests/benchmarks/old/dof_handler_timing.cc +++ /dev/null @@ -1,215 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2007 - 2013 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -// Refinement base: perform REF/dim refinements -// REF 21 needs 9GBytes on a 64 Bit architecture - -#ifdef DEBUG -#define REF 6 -#else -#define REF 21 -#endif - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace dealii; - -template -void indices (const DoFHandler &dof, unsigned int repeat) -{ - typedef typename DoFHandler::active_cell_iterator I; - - std::vector dofs(dof.get_fe().dofs_per_cell); - const I end = dof.end(); - - for (unsigned int k=0; kget_dof_indices(dofs); -} - - -template -void fevalues (const DoFHandler &dof, - UpdateFlags updates) -{ - typedef typename DoFHandler::active_cell_iterator I; - const I end = dof.end(); - - QGauss quadrature(5); - MappingQ1 mapping; - FEValues fe(mapping, dof.get_fe(), quadrature, updates); - - for (I i=dof.begin_active(); i!=end; ++i) - fe.reinit(i); -} - - -template -void fefacevalues (const DoFHandler &dof, - UpdateFlags updates) -{ - typedef typename DoFHandler::active_cell_iterator I; - const I end = dof.end(); - - QGauss quadrature(5); - MappingQ1 mapping; - FEFaceValues fe(mapping, dof.get_fe(), quadrature, updates); - - for (I i=dof.begin_active(); i!=end; ++i) - for (unsigned int f=0; f::faces_per_cell; ++f) - fe.reinit(i, f); -} - - -template -void check_mapping (const DoFHandler &dof) -{ - deallog.push("cell"); - fevalues(dof, update_q_points); - deallog << "qpoints" << std::endl; - fevalues(dof, update_JxW_values); - deallog << "JxW" << std::endl; - fevalues(dof, update_q_points | update_JxW_values); - deallog << "qpoints|JxW" << std::endl; - deallog.pop(); - - deallog.push("face"); - fefacevalues(dof, update_quadrature_points); - deallog << "qpoints" << std::endl; - fefacevalues(dof, update_JxW_values); - deallog << "JxW" << std::endl; - fefacevalues(dof, update_normal_vectors); - deallog << "normals" << std::endl; - fefacevalues(dof, update_q_points | update_JxW_values | update_normal_vectors); - deallog << "qpoints|JxW|normals" << std::endl; - deallog.pop(); - -} - - -template -void check_values (const DoFHandler &dof) -{ - indices(dof, 100); - deallog << "Index*100" << std::endl; - deallog.push("cell"); - fevalues(dof, update_values); - deallog << "values" << std::endl; - fevalues(dof, update_gradients); - deallog << "gradients" << std::endl; - fevalues(dof, update_values | update_JxW_values); - deallog << "values|JxW" << std::endl; - fevalues(dof, update_values | update_gradients - | update_q_points | update_JxW_values); - deallog << "values|gradients|qpoints|JxW" << std::endl; -// fevalues(dof, update_values | update_gradients | update_second_derivatives -// | update_q_points | update_JxW_values); -// deallog << "values|gradients|2nds|qpoints|JxW" << std::endl; - deallog.pop(); - - deallog.push("face"); - fefacevalues(dof, update_values); - deallog << "values" << std::endl; - fefacevalues(dof, update_gradients); - deallog << "gradients" << std::endl; - fefacevalues(dof, update_values | update_JxW_values); - deallog << "values|JxW" << std::endl; - fefacevalues(dof, update_values | update_gradients - | update_q_points | update_JxW_values | update_normal_vectors); - deallog << "values|gradients|qpoints|JxW|normals" << std::endl; -// fefacevalues(dof, update_values | update_gradients | update_second_derivatives -// | update_q_points | update_JxW_values); -// deallog << "values|gradients|2nds|qpoints|JxW" << std::endl; - deallog.pop(); -} - - -template -void check_q () -{ - deallog << "Mesh" << std::endl; - Triangulation tr; - GridGenerator::hyper_cube(tr); - tr.refine_global(REF/dim); - - for (unsigned int i=1; i<5; ++i) - { - FE_Q q(i); - deallog.push(q.get_name()); - deallog << "Dofs per cell " << q.dofs_per_cell << std::endl; - DoFHandler dof(tr); - dof.distribute_dofs(q); - deallog << "Dofs " << dof.n_dofs() << std::endl; - - if (i==1) - check_mapping(dof); - check_values(dof); - deallog.pop(); - } -} - - -template -void check_sys () -{ - deallog << "Mesh" << std::endl; - Triangulation tr; - GridGenerator::hyper_cube(tr); - tr.refine_global(REF/dim); - - FE_Q q(1); - - for (unsigned int i=1; i<5; ++i) - { - FESystem fe(q,1< dof(tr); - dof.distribute_dofs(fe); - deallog << "Dofs " << dof.n_dofs() << std::endl; - - if (i==1) - check_mapping(dof); - check_values(dof); - deallog.pop(); - } -} - - -int main() -{ - std::ofstream out("dof_handler_timing/output"); - deallog.attach(out); - deallog.log_execution_time(true); - deallog.log_time_differences(true); - check_q<2>(); - check_sys<2>(); - check_q<3>(); - check_sys<3>(); -} - diff --git a/tests/benchmarks/old/dof_handler_timing_01.cc b/tests/benchmarks/old/dof_handler_timing_01.cc deleted file mode 100644 index 2dfada080e..0000000000 --- a/tests/benchmarks/old/dof_handler_timing_01.cc +++ /dev/null @@ -1,116 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2007 - 2013 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace dealii; - - -// Fill dof handlers for different elements and see how large they get. -template -void check_dofs(DOF &dof, unsigned int repeat) -{ - FE_Q q1(1); - - deallog << "Start" << std::endl; - for (unsigned int i=1; i<6; ++i) - { - FE_Q fe(i); - for (unsigned int i=0; i fe(q1, 1 << i); - for (unsigned int i=0; i -void check (bool local) -{ - deallog << "Start mesh" << std::endl; - Triangulation tr(Triangulation::maximum_smoothing); - GridGenerator::hyper_cube(tr); - - if (local) - for (unsigned int i=0; i<1800/dim; ++i) - { - tr.last_active()->set_refine_flag(); - tr.execute_coarsening_and_refinement(); - } - else - tr.refine_global(21/dim); - - deallog << "Levels " << tr.n_levels() << " Cells "<< tr.n_cells() - << std::endl - << "active " << tr.n_active_cells() - << " memory " << tr.memory_consumption()/1.e6 - << std::endl; - - deallog.push("DoF"); - DoFHandler dof(tr); - check_dofs(dof, (local ? 100 : 10)); - deallog.pop(); - deallog.push("MGDoF"); - DoFHandler mgdof(tr); - check_dofs(mgdof, (local ? 100 : 10)); - deallog.pop(); -} - -int main() -{ - std::ofstream out("dof_handler_timing_01/output"); - deallog.log_execution_time(true); - deallog.log_time_differences(true); - deallog.attach(out); - deallog.push("local"); - deallog.push("2d"); - check<2>(true); - deallog.pop(); - deallog.push("3d"); - check<3>(true); - deallog.pop(); - deallog.pop(); - deallog.push("global"); - deallog.push("2d"); - check<2>(false); - deallog.pop(); - deallog.push("3d"); - check<3>(false); - deallog.pop(); - deallog.pop(); -} - diff --git a/tests/benchmarks/old/patches/library/2007-05-30 b/tests/benchmarks/old/patches/library/2007-05-30 deleted file mode 100644 index ac5c4fec5a..0000000000 --- a/tests/benchmarks/old/patches/library/2007-05-30 +++ /dev/null @@ -1,91 +0,0 @@ -Index: deal.II/source/multigrid/mg_dof_handler.cc -=================================================================== ---- deal.II/source/multigrid/mg_dof_handler.cc (revision 14714) -+++ deal.II/source/multigrid/mg_dof_handler.cc (working copy) -@@ -2,7 +2,7 @@ - // $Id$ - // Version: $Name$ - // --// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors -+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors - // - // This file is subject to QPL and may not be distributed - // without copyright and license information. Please refer -@@ -11,7 +11,7 @@ - // - //--------------------------------------------------------------------------- - -- -+#include - #include - #include - #include -@@ -33,6 +33,7 @@ - - - /* ------------------------ MGVertexDoFs ----------------------------------- */ -+//TODO: This seems horrible memory fragmentation! - - template - DoFHandler::MGVertexDoFs::MGVertexDoFs () -@@ -124,6 +125,25 @@ - } - - -+template -+unsigned int -+DoFHandler::memory_consumption() const -+{ -+ unsigned int mem = DoFHandler::memory_consumption(); -+ for (unsigned int l=0;lmemory_consumption(); -+ -+ mem += MemoryConsumption::memory_consumption(*mg_faces); -+ -+ for (unsigned int l=0;l -Index: deal.II/include/multigrid/mg_dof_handler.h -=================================================================== ---- deal.II/include/multigrid/mg_dof_handler.h (revision 14714) -+++ deal.II/include/multigrid/mg_dof_handler.h (working copy) -@@ -2,7 +2,7 @@ - // $Id$ - // Version: $Name$ - // --// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors -+// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors - // - // This file is subject to QPL and may not be distributed - // without copyright and license information. Please refer -@@ -711,6 +711,21 @@ - */ - unsigned int n_dofs () const; - -+ /** -+ * Determine an estimate for the -+ * memory consumption (in bytes) -+ * of this object. -+ * -+ * This function is made virtual, -+ * since a dof handler object -+ * might be accessed through a -+ * pointers to thisr base class, -+ * although the actual object -+ * might be a derived class. -+ */ -+ virtual unsigned int memory_consumption () const; -+ -+ - /** - * Exception. - */ diff --git a/tests/benchmarks/old/performance.cc b/tests/benchmarks/old/performance.cc deleted file mode 100644 index b16d91d6cf..0000000000 --- a/tests/benchmarks/old/performance.cc +++ /dev/null @@ -1,141 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2013 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -#include "../tests.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - - -template -void performance (Triangulation &tr, - const Mapping &mapping, - const FiniteElement &fe, - const Quadrature &quadrature, - UpdateFlags flags) -{ - deallog << "Create dofs" << std::endl; - DoFHandler dof (tr); - dof.distribute_dofs (fe); - - deallog << "Create FEValues" << std::endl; - - FEValues val (mapping, fe, quadrature, flags); - - DoFHandler::active_cell_iterator cell = dof.begin_active(); - DoFHandler::active_cell_iterator end = dof.end(); - - deallog << "Loop" << std::endl; - - for (; cell != end ; ++cell) - val.reinit(cell); - - deallog << "End" << std::endl; -} - -template -void loop (std::vector *> elements, - const Mapping &mapping, - Triangulation &tr) -{ - QGauss gauss (4); - - typename std::vector *>::iterator elementp = elements.begin (); - typename std::vector *>::iterator end = elements.end (); - - for (; elementp != end; ++elementp) - { - const FiniteElement &element = **elementp; - - char dofs[20]; - std::ostrstream ost (dofs, 19); - ost << element.n_dofs_per_cell() << std::ends; - - deallog.push(dofs); - - deallog.push("points"); - performance (tr, mapping, element, gauss, update_q_points); - deallog.pop(); - - deallog.push("values"); - performance (tr, mapping, element, gauss, update_values); - deallog.pop(); - - deallog.push("grads-"); - performance (tr, mapping, element, gauss, update_gradients); - deallog.pop(); - - deallog.push("2nd---"); - performance (tr, mapping, element, gauss, update_second_derivatives); - deallog.pop(); - - deallog.push("matrix"); - performance (tr, mapping, element, gauss, update_q_points - | update_JxW_values - | update_values - | update_gradients); - deallog.pop(); - - deallog.push("all---"); - performance (tr, mapping, element, gauss, update_q_points - | update_JxW_values - | update_values - | update_gradients - | update_second_derivatives); - deallog.pop(); - deallog.pop(); - } -} - -int main () -{ - std::ofstream of ("performance.log"); - deallog.attach (of); - deallog.log_execution_time(true); - deallog.log_time_differences(true); - Triangulation<2> tr; - GridGenerator::hyper_ball (tr); - tr.refine_global (8); - - MappingCartesian<2> cartesian; - MappingQ1<2> mapping; - MappingQ<2> mappingq1(1); - MappingQ<2> mappingq2(2); - std::vector*> el2d; - el2d.push_back (new FE_Q<2> (1)); - el2d.push_back (new FE_Q<2> (2)); - el2d.push_back (new FE_Q<2> (3)); - el2d.push_back (new FE_Q<2> (4)); - - loop (el2d, cartesian, tr); - loop (el2d, mapping, tr); - loop (el2d, mappingq1, tr); - loop (el2d, mappingq2, tr); -} diff --git a/tests/benchmarks/old/performance.pl b/tests/benchmarks/old/performance.pl deleted file mode 100644 index ed3ffff689..0000000000 --- a/tests/benchmarks/old/performance.pl +++ /dev/null @@ -1,27 +0,0 @@ -# -# Copyright (C) 2001 by the deal.II authors -# -# -# Make gnuplot table from output of performance -# -# Usage: -# perl performance.pl performance.log -# -# order of values per line is -# 1. dofs per cell -# 2. Time for generating quadrature points -# 3. Time for generating shape values -# 4. Time for generating shape gradients -# 5. Time for generating shape second derivatives -# 6. Time for generating typical mix for matrix generation -# 7. Time for generating all of the above - -while (<>) -{ - if (/(\d+\.\d+):DEAL:(\d+):([^:]+)::End/) - { - printf "%3d", $2 if ($3 eq "points"); - printf "\t%f", $1; - print "\n" if ($3 eq "all---"); - } -} diff --git a/tests/benchmarks/old/spec2006-447.dealII.cc b/tests/benchmarks/old/spec2006-447.dealII.cc deleted file mode 100644 index f009095a5b..0000000000 --- a/tests/benchmarks/old/spec2006-447.dealII.cc +++ /dev/null @@ -1,4279 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2002 - 2014 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - - -// Start out with well known things... - - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#ifdef HAVE_STD_STRINGSTREAM -# include -#else -# include -#endif - -using namespace dealii; - - -/** - * The number of iterations the program shall do. This is given as a - * parameter to the executable. - */ -unsigned int n_steps = 0; -/** - * The present step. - */ -unsigned int step; - -/** - * Declare the coefficient in front - * of the Laplace operator of the - * equation. - */ -template -class LaplaceCoefficient : public Function -{ -public: - virtual double value (const Point &p, - const unsigned int) const - { - return 1.+p*p; - } -}; - - -/** - * Same for the coefficient in front - * of the mass term. - */ -template -class MassCoefficient : public Function -{ -public: - virtual double value (const Point &p, - const unsigned int) const - { - double factor = 1; - for (unsigned int d=0; d - class EvaluationBase - { - public: - virtual ~EvaluationBase (); - - void set_refinement_cycle (const unsigned int refinement_cycle); - - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const = 0; - - unsigned int refinement_cycle; - }; - - - template - EvaluationBase::~EvaluationBase () - {} - - - - template - void - EvaluationBase::set_refinement_cycle (const unsigned int step) - { - refinement_cycle = step; - } - - - // @sect4{The PointValueEvaluation class} - template - class PointValueEvaluation : public EvaluationBase - { - public: - PointValueEvaluation (const Point &evaluation_point); - - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; - - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); - - const Point evaluation_point; - }; - - - template - PointValueEvaluation:: - PointValueEvaluation (const Point &evaluation_point) - : - evaluation_point (evaluation_point) - {} - - - - template - void - PointValueEvaluation:: - operator () (const DoFHandler &dof_handler, - const Vector &solution) const - { - double point_value = 1e20; - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - bool evaluation_point_found = false; - for (; (cell!=endc) && !evaluation_point_found; ++cell) - for (unsigned int vertex=0; - vertex::vertices_per_cell; - ++vertex) - if (cell->vertex(vertex).distance (evaluation_point) - < - cell->diameter() * 1e-8) - { - point_value = solution(cell->vertex_dof_index(vertex,0)); - - evaluation_point_found = true; - break; - }; - - AssertThrow (evaluation_point_found, - ExcEvaluationPointNotFound(evaluation_point)); - - std::cout << " Point value=" << point_value - << std::endl; - } - - - // @sect4{The PointXDerivativeEvaluation class} - - // Besides the class implementing - // the evaluation of the solution - // at one point, we here provide - // one which evaluates the gradient - // at a grid point. Since in - // general the gradient of a finite - // element function is not - // continuous at a vertex, we have - // to be a little bit more careful - // here. What we do is to loop over - // all cells, even if we have found - // the point already on one cell, - // and use the mean value of the - // gradient at the vertex taken - // from all adjacent cells. - // - // Given the interface of the - // ``PointValueEvaluation'' class, - // the declaration of this class - // provides little surprise, and - // neither does the constructor: - template - class PointXDerivativeEvaluation : public EvaluationBase - { - public: - PointXDerivativeEvaluation (const Point &evaluation_point); - - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; - - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); - - const Point evaluation_point; - }; - - - template - PointXDerivativeEvaluation:: - PointXDerivativeEvaluation (const Point &evaluation_point) - : - evaluation_point (evaluation_point) - {} - - - // The more interesting things - // happen inside the function doing - // the actual evaluation: - template - void - PointXDerivativeEvaluation:: - operator () (const DoFHandler &dof_handler, - const Vector &solution) const - { - // This time initialize the - // return value with something - // useful, since we will have to - // add up a number of - // contributions and take the - // mean value afterwards... - double point_derivative = 0; - - // ...then have some objects of - // which the meaning wil become - // clear below... - QTrapez vertex_quadrature; - MappingQ mapping (4); - FEValues fe_values (mapping, dof_handler.get_fe(), - vertex_quadrature, - update_gradients | update_q_points); - std::vector > - solution_gradients (vertex_quadrature.n_quadrature_points); - - // ...and next loop over all cells - // and their vertices, and count - // how often the vertex has been - // found: - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - unsigned int evaluation_point_hits = 0; - for (; cell!=endc; ++cell) - for (unsigned int vertex=0; - vertex::vertices_per_cell; - ++vertex) - if (cell->vertex(vertex) == evaluation_point) - { - // Things are now no more - // as simple, since we - // can't get the gradient - // of the finite element - // field as before, where - // we simply had to pick - // one degree of freedom - // at a vertex. - // - // Rather, we have to - // evaluate the finite - // element field on this - // cell, and at a certain - // point. As you know, - // evaluating finite - // element fields at - // certain points is done - // through the - // ``FEValues'' class, so - // we use that. The - // question is: the - // ``FEValues'' object - // needs to be a given a - // quadrature formula and - // can then compute the - // values of finite - // element quantities at - // the quadrature - // points. Here, we don't - // want to do quadrature, - // we simply want to - // specify some points! - // - // Nevertheless, the same - // way is chosen: use a - // special quadrature - // rule with points at - // the vertices, since - // these are what we are - // interested in. The - // appropriate rule is - // the trapezoidal rule, - // so that is the reason - // why we used that one - // above. - // - // Thus: initialize the - // ``FEValues'' object on - // this cell, - fe_values.reinit (cell); - // and extract the - // gradients of the - // solution vector at the - // vertices: - fe_values.get_function_grads (solution, - solution_gradients); - - // Now we have the - // gradients at all - // vertices, so pick out - // that one which belongs - // to the evaluation - // point (note that the - // order of vertices is - // not necessarily the - // same as that of the - // quadrature points): - unsigned int q_point = 0; - for (; q_point 0, - ExcEvaluationPointNotFound(evaluation_point)); - - // We have simply summed up the - // contributions of all adjacent - // cells, so we still have to - // compute the mean value. Once - // this is done, report the status: - point_derivative /= evaluation_point_hits; - std::cout << " Point x-derivative=" << point_derivative - << std::endl; - } - - - - // @sect4{The GridOutput class} - - // Since this program has a more - // difficult structure (it computed - // a dual solution in addition to a - // primal one), writing out the - // solution is no more done by an - // evaluation object since we want - // to write both solutions at once - // into one file, and that requires - // some more information than - // available to the evaluation - // classes. - // - // However, we also want to look at - // the grids generated. This again - // can be done with one such - // class. Its structure is analog - // to the ``SolutionOutput'' class - // of the previous example program, - // so we do not discuss it here in - // more detail. Furthermore, - // everything that is used here has - // already been used in previous - // example programs. - template - class GridOutput : public EvaluationBase - { - public: - GridOutput (const std::string &output_name_base); - - virtual void operator () (const DoFHandler &dof_handler, - const Vector &solution) const; - - const std::string output_name_base; - }; - - - template - GridOutput:: - GridOutput (const std::string &output_name_base) - : - output_name_base (output_name_base) - {} - - - template - void - GridOutput::operator () (const DoFHandler &dof_handler, - const Vector &/*solution*/) const - { -#ifdef HAVE_STD_STRINGSTREAM - std::ostringstream filename; -#else - std::ostrstream filename; -#endif - filename << "spec2006-447.dealII/" - << output_name_base << "-" - << this->refinement_cycle - << ".eps" - << std::ends; -#ifdef HAVE_STD_STRINGSTREAM - std::ofstream out (filename.str().c_str()); -#else - std::ofstream out (filename.str()); -#endif - - GridOut().write_eps (dof_handler.get_tria(), out); - } -} - - -// @sect3{The Laplace solver classes} - -// Next are the actual solver -// classes. Again, we discuss only -// the differences to the previous -// program. -namespace LaplaceSolver -{ - // Before everything else, - // forward-declare one class that - // we will have later, since we - // will want to make it a friend of - // some of the classes that follow, - // which requires the class to be - // known: - template class WeightedResidual; - - - // @sect4{The Laplace solver base class} - - // This class is almost unchanged, - // with the exception that it - // declares two more functions: - // ``output_solution'' will be used - // to generate output files from - // the actual solutions computed by - // derived classes, and the - // ``set_refinement_cycle'' - // function by which the testing - // framework sets the number of the - // refinement cycle to a local - // variable in this class; this - // number is later used to generate - // filenames for the solution - // output. - template - class Base - { - public: - Base (Triangulation &coarse_grid); - virtual ~Base (); - - virtual void solve_problem () = 0; - virtual void postprocess (const Evaluation::EvaluationBase &postprocessor) const = 0; - virtual void refine_grid () = 0; - virtual unsigned int n_dofs () const = 0; - - virtual void set_refinement_cycle (const unsigned int cycle); - - virtual void output_solution () const = 0; - - - const SmartPointer > triangulation; - - unsigned int refinement_cycle; - }; - - - template - Base::Base (Triangulation &coarse_grid) - : - triangulation (&coarse_grid) - {} - - - template - Base::~Base () - {} - - - - template - void - Base::set_refinement_cycle (const unsigned int cycle) - { - refinement_cycle = cycle; - } - - - // @sect4{The Laplace Solver class} - - // Likewise, the ``Solver'' class - // is entirely unchanged and will - // thus not be discussed. - template - class Solver : public virtual Base - { - public: - Solver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &boundary_values); - virtual - ~Solver (); - - virtual - void - solve_problem (); - - virtual - void - postprocess (const Evaluation::EvaluationBase &postprocessor) const; - - virtual - unsigned int - n_dofs () const; - - - const SmartPointer > fe; - const SmartPointer > quadrature; - const SmartPointer > face_quadrature; - DoFHandler dof_handler; - Vector solution; - const SmartPointer > boundary_values; - - virtual void assemble_rhs (Vector &rhs) const = 0; - - - struct LinearSystem - { - LinearSystem (const DoFHandler &dof_handler); - - void solve (Vector &solution) const; - - ConstraintMatrix hanging_node_constraints; - SparsityPattern sparsity_pattern; - SparseMatrix matrix; - Vector rhs; - }; - - void - assemble_linear_system (LinearSystem &linear_system); - - void - assemble_matrix (LinearSystem &linear_system, - const typename DoFHandler::active_cell_iterator &begin_cell, - const typename DoFHandler::active_cell_iterator &end_cell, - Threads::Mutex &mutex) const; - }; - - - - template - Solver::Solver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &boundary_values) - : - Base (triangulation), - fe (&fe), - quadrature (&quadrature), - face_quadrature (&face_quadrature), - dof_handler (triangulation), - boundary_values (&boundary_values) - {} - - - template - Solver::~Solver () - { - dof_handler.clear (); - } - - - template - void - Solver::solve_problem () - { - dof_handler.distribute_dofs (*fe); - solution.reinit (dof_handler.n_dofs()); - - LinearSystem linear_system (dof_handler); - assemble_linear_system (linear_system); - linear_system.solve (solution); - } - - - template - void - Solver:: - postprocess (const Evaluation::EvaluationBase &postprocessor) const - { - postprocessor (dof_handler, solution); - } - - - template - unsigned int - Solver::n_dofs () const - { - return dof_handler.n_dofs(); - } - - - template - void - Solver::assemble_linear_system (LinearSystem &linear_system) - { - typedef - typename DoFHandler::active_cell_iterator - active_cell_iterator; - - const unsigned int n_threads = multithread_info.n_default_threads; - std::vector > - thread_ranges - = Threads::split_range (dof_handler.begin_active (), - dof_handler.end (), - n_threads); - - Threads::Mutex mutex; - Threads::ThreadGroup<> threads; - for (unsigned int thread=0; thread::assemble_matrix, *this, - linear_system, - thread_ranges[thread].first, - thread_ranges[thread].second, - mutex); - - assemble_rhs (linear_system.rhs); - linear_system.hanging_node_constraints.condense (linear_system.rhs); - - std::map boundary_value_map; - VectorTools::interpolate_boundary_values (dof_handler, - 0, - *boundary_values, - boundary_value_map); - - threads.join_all (); - linear_system.hanging_node_constraints.condense (linear_system.matrix); - - MatrixTools::apply_boundary_values (boundary_value_map, - linear_system.matrix, - solution, - linear_system.rhs); - } - - - template - void - Solver::assemble_matrix (LinearSystem &linear_system, - const typename DoFHandler::active_cell_iterator &begin_cell, - const typename DoFHandler::active_cell_iterator &end_cell, - Threads::Mutex &mutex) const - { - MappingQ mapping (4); - FEValues fe_values (mapping, *fe, *quadrature, - UpdateFlags(update_gradients | update_values | - update_q_points | - update_JxW_values)); - - const unsigned int dofs_per_cell = fe->dofs_per_cell; - const unsigned int n_q_points = quadrature->n_quadrature_points; - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - - std::vector local_dof_indices (dofs_per_cell); - - std::vector laplace_coefficients (fe_values.n_quadrature_points); - std::vector mass_coefficients (fe_values.n_quadrature_points); - - - for (typename DoFHandler::active_cell_iterator cell=begin_cell; - cell!=end_cell; ++cell) - { - cell_matrix = 0; - - fe_values.reinit (cell); - - LaplaceCoefficient().value_list (fe_values.get_quadrature_points(), - laplace_coefficients); - MassCoefficient().value_list (fe_values.get_quadrature_points(), - mass_coefficients); - - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - Threads::Mutex::ScopedLock lock (mutex); - for (unsigned int i=0; i - Solver::LinearSystem:: - LinearSystem (const DoFHandler &dof_handler) - { - hanging_node_constraints.clear (); - - void (*mhnc_p) (const DoFHandler &, - ConstraintMatrix &) - = &DoFTools::make_hanging_node_constraints; - - Threads::Thread<> - mhnc_thread = Threads::new_thread (mhnc_p,dof_handler, hanging_node_constraints); - - // make sparsity pattern. since - // in 3d the usual way just blows - // the roof w.r.t. memory - // consumption, use the detour - // via a compressed sparsity - // pattern that we later copy - // over - CompressedSparsityPattern csp(dof_handler.n_dofs(), - dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, csp); - - mhnc_thread.join (); - hanging_node_constraints.close (); - hanging_node_constraints.condense (csp); - - sparsity_pattern.copy_from(csp); - matrix.reinit (sparsity_pattern); - rhs.reinit (dof_handler.n_dofs()); - } - - - - template - void - Solver::LinearSystem::solve (Vector &solution) const - { - SolverControl solver_control (solution.size(), 1e-6); - PrimitiveVectorMemory<> vector_memory; - SolverCG<> cg (solver_control, vector_memory); - - PreconditionJacobi<> preconditioner; - preconditioner.initialize(matrix); - - cg.solve (matrix, solution, rhs, preconditioner); - - hanging_node_constraints.distribute (solution); - } - - - - - // @sect4{The PrimalSolver class} - - // The ``PrimalSolver'' class is - // also mostly unchanged except for - // overloading the functions - // ``solve_problem'', ``n_dofs'', - // and ``postprocess'' of the base - // class, and implementing the - // ``output_solution'' - // function. These overloaded - // functions do nothing particular - // besides calling the functions of - // the base class -- that seems - // superfluous, but works around a - // bug in a popular compiler which - // requires us to write such - // functions for the following - // scenario: Besides the - // ``PrimalSolver'' class, we will - // have a ``DualSolver'', both - // derived from ``Solver''. We will - // then have a final classes which - // derived from these two, which - // will then have two instances of - // the ``Solver'' class as its base - // classes. If we want, for - // example, the number of degrees - // of freedom of the primal solver, - // we would have to indicate this - // like so: - // ``PrimalSolver::n_dofs()''. - // However, the compiler does not - // accept this since the ``n_dofs'' - // function is actually from a base - // class of the ``PrimalSolver'' - // class, so we have to inject the - // name from the base to the - // derived class using these - // additional functions. - // - // Regarding the implementation of - // the ``output_solution'' - // function, we keep the - // ``GlobalRefinement'' and - // ``RefinementKelly'' classes in - // this program, and they can then - // rely on the default - // implementation of this function - // which simply outputs the primal - // solution. The class implementing - // dual weighted error estimators - // will overload this function - // itself, to also output the dual - // solution. - // - // Except for this, the class is - // unchanged with respect to the - // previous example. - template - class PrimalSolver : public Solver - { - public: - PrimalSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values); - - virtual - void solve_problem (); - - virtual - unsigned int n_dofs () const; - - virtual - void postprocess (const Evaluation::EvaluationBase &postprocessor) const; - - virtual - void output_solution () const; - - - const SmartPointer > rhs_function; - virtual void assemble_rhs (Vector &rhs) const; - - // Now, in order to work around - // some problems in one of the - // compilers this library can - // be compiled with, we will - // have to use some - // workarounds. This will - // require that we declare a - // class that is actually - // derived from the present - // one, as a friend (strange as - // that seems). The full - // rationale will be explained - // below. - friend class WeightedResidual; - }; - - - template - PrimalSolver:: - PrimalSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values) - : - Base (triangulation), - Solver (triangulation, fe, - quadrature, face_quadrature, - boundary_values), - rhs_function (&rhs_function) - {} - - - template - void - PrimalSolver::solve_problem () - { - Solver::solve_problem (); - } - - - - template - unsigned int - PrimalSolver::n_dofs() const - { - return Solver::n_dofs(); - } - - - template - void - PrimalSolver:: - postprocess (const Evaluation::EvaluationBase &postprocessor) const - { - Solver::postprocess(postprocessor); - } - - - template - void - PrimalSolver::output_solution () const - { - abort (); - } - - - - template - void - PrimalSolver:: - assemble_rhs (Vector &rhs) const - { - MappingQ mapping (4); - FEValues fe_values (mapping, *this->fe, *this->quadrature, - UpdateFlags(update_values | - update_q_points | - update_JxW_values)); - - const unsigned int dofs_per_cell = this->fe->dofs_per_cell; - const unsigned int n_q_points = this->quadrature->n_quadrature_points; - - Vector cell_rhs (dofs_per_cell); - std::vector rhs_values (n_q_points); - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell = this->dof_handler.begin_active(), - endc = this->dof_handler.end(); - for (; cell!=endc; ++cell) - { - cell_rhs = 0; - - fe_values.reinit (cell); - - rhs_function->value_list (fe_values.get_quadrature_points(), - rhs_values); - - for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); - for (unsigned int i=0; i - class RefinementGlobal : public PrimalSolver - { - public: - RefinementGlobal (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values); - - virtual void refine_grid (); - }; - - - - template - RefinementGlobal:: - RefinementGlobal (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values) - : - Base (coarse_grid), - PrimalSolver (coarse_grid, fe, quadrature, - face_quadrature, rhs_function, - boundary_values) - {} - - - - template - void - RefinementGlobal::refine_grid () - { - this->triangulation->refine_global (1); - } - - - - template - class RefinementKelly : public PrimalSolver - { - public: - RefinementKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values); - - virtual void refine_grid (); - }; - - - - template - RefinementKelly:: - RefinementKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values) - : - Base (coarse_grid), - PrimalSolver (coarse_grid, fe, quadrature, - face_quadrature, - rhs_function, boundary_values) - {} - - - - template - void - RefinementKelly::refine_grid () - { - Vector estimated_error_per_cell (this->triangulation->n_active_cells()); - KellyErrorEstimator::estimate (this->dof_handler, - QGauss(3), - typename FunctionMap::type(), - this->solution, - estimated_error_per_cell); - GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation, - estimated_error_per_cell, - 0.2, 0.02); - this->triangulation->execute_coarsening_and_refinement (); - } - - - - // @sect4{The RefinementWeightedKelly class} - - // This class is a variant of the - // previous one, in that it allows - // to weight the refinement - // indicators we get from the - // library's Kelly indicator by - // some function. We include this - // class since the goal of this - // example program is to - // demonstrate automatic refinement - // criteria even for complex output - // quantities such as point values - // or stresses. If we did not solve - // a dual problem and compute the - // weights thereof, we would - // probably be tempted to give a - // hand-crafted weighting to the - // indicators to account for the - // fact that we are going to - // evaluate these quantities. This - // class accepts such a weighting - // function as argument to its - // constructor: - template - class RefinementWeightedKelly : public PrimalSolver - { - public: - RefinementWeightedKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values, - const Function &weighting_function); - - virtual void refine_grid (); - - - const SmartPointer > weighting_function; - }; - - - - template - RefinementWeightedKelly:: - RefinementWeightedKelly (Triangulation &coarse_grid, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values, - const Function &weighting_function) - : - Base (coarse_grid), - PrimalSolver (coarse_grid, fe, quadrature, - face_quadrature, - rhs_function, boundary_values), - weighting_function (&weighting_function) - {} - - - - // Now, here comes the main - // function, including the - // weighting: - template - void - RefinementWeightedKelly::refine_grid () - { - // First compute some residual - // based error indicators for all - // cells by a method already - // implemented in the - // library. What exactly is - // computed can be read in the - // documentation of that class. - Vector estimated_error (this->triangulation->n_active_cells()); - KellyErrorEstimator::estimate (this->dof_handler, - *this->face_quadrature, - typename FunctionMap::type(), - this->solution, - estimated_error); - - // Now we are going to weight - // these indicators by the value - // of the function given to the - // constructor: - typename DoFHandler::active_cell_iterator - cell = this->dof_handler.begin_active(), - endc = this->dof_handler.end(); - for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index) - estimated_error(cell_index) - *= weighting_function->value (cell->center()); - - GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation, - estimated_error, - 0.2, 0.02); - this->triangulation->execute_coarsening_and_refinement (); - } - -} - - -// @sect3{Equation data} -// -// In this example program, we work -// with the same data sets as in the -// previous one, but as it may so -// happen that someone wants to run -// the program with different -// boundary values and right hand side -// functions, or on a different grid, -// we show a simple technique to do -// exactly that. For more clarity, we -// furthermore pack everything that -// has to do with equation data into -// a namespace of its own. -// -// The underlying assumption is that -// this is a research program, and -// that there we often have a number -// of test cases that consist of a -// domain, a right hand side, -// boundary values, possibly a -// specified coefficient, and a -// number of other parameters. They -// often vary all at the same time -// when shifting from one example to -// another. To make handling such -// sets of problem description -// parameters simple is the goal of -// the following. -// -// Basically, the idea is this: let -// us have a structure for each set -// of data, in which we pack -// everything that describes a test -// case: here, these are two -// subclasses, one called -// ``BoundaryValues'' for the -// boundary values of the exact -// solution, and one called -// ``RightHandSide'', and then a way -// to generate the coarse grid. Since -// the solution of the previous -// example program looked like curved -// ridges, we use this name here for -// the enclosing class. Note that the -// names of the two inner classes -// have to be the same for all -// enclosing test case classes, and -// also that we have attached the -// dimension template argument to the -// enclosing class rather than to the -// inner ones, to make further -// processing simpler. (From a -// language viewpoint, a namespace -// would be better to encapsulate -// these inner classes, rather than a -// structure. However, namespaces -// cannot be given as template -// arguments, so we use a structure -// to allow a second object to select -// from within its given -// argument. The enclosing structure, -// of course, has no member variables -// apart from the classes it -// declares, and a static function to -// generate the coarse mesh; it will -// in general never be instantiated.) -// -// The idea is then the following -// (this is the right time to also -// take a brief look at the code -// below): we can generate objects -// for boundary values and -// right hand side by simply giving -// the name of the outer class as a -// template argument to a class which -// we call here ``Data::SetUp'', and -// it then creates objects for the -// inner classes. In this case, to -// get all that characterizes the -// curved ridge solution, we would -// simply generate an instance of -// ``Data::SetUp'', -// and everything we need to know -// about the solution would be static -// member variables and functions of -// that object. -// -// This approach might seem like -// overkill in this case, but will -// become very handy once a certain -// set up is not only characterized -// by Dirichlet boundary values and a -// right hand side function, but in -// addition by material properties, -// Neumann values, different boundary -// descriptors, etc. In that case, -// the ``SetUp'' class might consist -// of a dozen or more objects, and -// each descriptor class (like the -// ``CurvedRidges'' class below) -// would have to provide them. Then, -// you will be happy to be able to -// change from one set of data to -// another by only changing the -// template argument to the ``SetUp'' -// class at one place, rather than at -// many. -// -// With this framework for different -// test cases, we are almost -// finished, but one thing remains: -// by now we can select statically, -// by changing one template argument, -// which data set to choose. In order -// to be able to do that dynamically, -// i.e. at run time, we need a base -// class. This we provide in the -// obvious way, see below, with -// virtual abstract functions. It -// forces us to introduce a second -// template parameter ``dim'' which -// we need for the base class (which -// could be avoided using some -// template magic, but we omit that), -// but that's all. -// -// Adding new testcases is now -// simple, you don't have to touch -// the framework classes, only a -// structure like the -// ``CurvedRidges'' one is needed. -namespace Data -{ - // @sect4{The SetUpBase and SetUp classes} - - // Based on the above description, - // the ``SetUpBase'' class then - // looks as follows. To allow using - // the ``SmartPointer'' class with - // this class, we derived from the - // ``Subscriptor'' class. - template - struct SetUpBase : public Subscriptor - { - virtual - const Function &get_boundary_values () const = 0; - - virtual - const Function &get_right_hand_side () const = 0; - - virtual - void create_coarse_grid (Triangulation &coarse_grid) const = 0; - }; - - - // And now for the derived class - // that takes the template argument - // as explained above. For some - // reason, C++ requires us to - // define a constructor (which - // maybe empty), as otherwise a - // warning is generated that some - // data is not initialized. - // - // Here we pack the data elements - // into private variables, and - // allow access to them through the - // methods of the base class. - template - struct SetUp : public SetUpBase - { - SetUp () {}; - - virtual - const Function &get_boundary_values () const; - - virtual - const Function &get_right_hand_side () const; - - - virtual - void create_coarse_grid (Triangulation &coarse_grid) const; - - - static const typename Traits::BoundaryValues boundary_values; - static const typename Traits::RightHandSide right_hand_side; - }; - - // We have to provide definitions - // for the static member variables - // of the above class: - template - const typename Traits::BoundaryValues SetUp::boundary_values; - template - const typename Traits::RightHandSide SetUp::right_hand_side; - - // And definitions of the member - // functions: - template - const Function & - SetUp::get_boundary_values () const - { - return boundary_values; - } - - - template - const Function & - SetUp::get_right_hand_side () const - { - return right_hand_side; - } - - - template - void - SetUp:: - create_coarse_grid (Triangulation &coarse_grid) const - { - Traits::create_coarse_grid (coarse_grid); - } - - - // @sect4{The CurvedRidges class} - - // The class that is used to - // describe the boundary values and - // right hand side of the ``curved - // ridge'' problem already used in - // the step-13 example program is - // then like so: - template - struct CurvedRidges - { - class BoundaryValues : public Function - { - public: - BoundaryValues () : Function () {}; - - virtual double value (const Point &p, - const unsigned int component) const; - }; - - - class RightHandSide : public Function - { - public: - RightHandSide () : Function () {}; - - virtual double value (const Point &p, - const unsigned int component) const; - }; - - static - void - create_coarse_grid (Triangulation &coarse_grid); - }; - - - template - double - CurvedRidges::BoundaryValues:: - value (const Point &p, - const unsigned int /*component*/) const - { - double q = p(0); - for (unsigned int i=1; i - double - CurvedRidges::RightHandSide::value (const Point &p, - const unsigned int /*component*/) const - { - double q = p(0); - for (unsigned int i=1; i - void - CurvedRidges:: - create_coarse_grid (Triangulation &coarse_grid) - { - GridGenerator::hyper_cube (coarse_grid, -1, 1); - coarse_grid.refine_global (2); - } - - - // @sect4{The Exercise_2_3 class} - - // This example program was written - // while giving practical courses - // for a lecture on adaptive finite - // element methods and duality - // based error estimates. For these - // courses, we had one exercise, - // which required to solve the - // Laplace equation with constant - // right hand side on a square - // domain with a square hole in the - // center, and zero boundary - // values. Since the implementation - // of the properties of this - // problem is so particularly - // simple here, lets do it. As the - // number of the exercise was 2.3, - // we take the liberty to retain - // this name for the class as well. - template - struct Exercise_2_3 - { - // We need a class to denote - // the boundary values of the - // problem. In this case, this - // is simple: it's the zero - // function, so don't even - // declare a class, just a - // typedef: - typedef ZeroFunction BoundaryValues; - - // Second, a class that denotes - // the right hand side. Since - // they are constant, just - // subclass the corresponding - // class of the library and be - // done: - class RightHandSide : public ConstantFunction - { - public: - RightHandSide () : ConstantFunction (1.) {}; - }; - - // Finally a function to - // generate the coarse - // grid. This is somewhat more - // complicated here, see - // immediately below. - static - void - create_coarse_grid (Triangulation &coarse_grid); - }; - - - // As stated above, the grid for - // this example is the square - // [-1,1]^2 with the square - // [-1/2,1/2]^2 as hole in it. We - // create the coarse grid as 4 - // times 4 cells with the middle - // four ones missing. - // - // Of course, the example has an - // extension to 3d, but since this - // function cannot be written in a - // dimension independent way we - // choose not to implement this - // here, but rather only specialize - // the template for dim=2. If you - // compile the program for 3d, - // you'll get a message from the - // linker that this function is not - // implemented for 3d, and needs to - // be provided. - // - // For the creation of this - // geometry, the library has no - // predefined method. In this case, - // the geometry is still simple - // enough to do the creation by - // hand, rather than using a mesh - // generator. - /* - template <> - void - Exercise_2_3<2>:: - create_coarse_grid (Triangulation<2> &coarse_grid) - { - // First define the space - // dimension, to allow those - // parts of the function that are - // actually dimension independent - // to use this variable. That - // makes it simpler if you later - // takes this as a starting point - // to implement the 3d version. - const unsigned int dim = 2; - - // Then have a list of - // vertices. Here, they are 24 (5 - // times 5, with the middle one - // omitted). It is probably best - // to draw a sketch here. Note - // that we leave the number of - // vertices open at first, but - // then let the compiler compute - // this number afterwards. This - // reduces the possibility of - // having the dimension to large - // and leaving the last ones - // uninitialized. - static const Point<2> vertices_1[] - = { Point<2> (-1., -1.), - Point<2> (-1./2, -1.), - Point<2> (0., -1.), - Point<2> (+1./2, -1.), - Point<2> (+1, -1.), - - Point<2> (-1., -1./2.), - Point<2> (-1./2, -1./2.), - Point<2> (0., -1./2.), - Point<2> (+1./2, -1./2.), - Point<2> (+1, -1./2.), - - Point<2> (-1., 0.), - Point<2> (-1./2, 0.), - Point<2> (+1./2, 0.), - Point<2> (+1, 0.), - - Point<2> (-1., 1./2.), - Point<2> (-1./2, 1./2.), - Point<2> (0., 1./2.), - Point<2> (+1./2, 1./2.), - Point<2> (+1, 1./2.), - - Point<2> (-1., 1.), - Point<2> (-1./2, 1.), - Point<2> (0., 1.), - Point<2> (+1./2, 1.), - Point<2> (+1, 1.) }; - const unsigned int - n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]); - - // From this static list of - // vertices, we generate an STL - // vector of the vertices, as - // this is the data type the - // library wants to see. - const std::vector > vertices (&vertices_1[0], - &vertices_1[n_vertices]); - - // Next, we have to define the - // cells and the vertices they - // contain. Here, we have 8 - // vertices, but leave the number - // open and let it be computed - // afterwards: - static const int cell_vertices[][GeometryInfo::vertices_per_cell] - = {{0, 1, 6,5}, - {1, 2, 7, 6}, - {2, 3, 8, 7}, - {3, 4, 9, 8}, - {5, 6, 11, 10}, - {8, 9, 13, 12}, - {10, 11, 15, 14}, - {12, 13, 18, 17}, - {14, 15, 20, 19}, - {15, 16, 21, 20}, - {16, 17, 22, 21}, - {17, 18, 23, 22}}; - const unsigned int - n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]); - - // Again, we generate a C++ - // vector type from this, but - // this time by looping over the - // cells (yes, this is - // boring). Additionally, we set - // the material indicator to zero - // for all the cells: - std::vector > cells (n_cells, CellData()); - for (unsigned int i=0; i::vertices_per_cell; - ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - }; - - // Finally pass all this - // information to the library to - // generate a triangulation. The - // last parameter may be used to - // pass information about - // non-zero boundary indicators - // at certain faces of the - // triangulation to the library, - // but we don't want that here, - // so we give an empty object: - coarse_grid.create_triangulation (vertices, - cells, - SubCellData()); - - // And since we want that the - // evaluation point (3/4,3/4) in - // this example is a grid point, - // we refine once globally: - coarse_grid.refine_global (1); - } - */ - - - template <> - void - Exercise_2_3<3>:: - create_coarse_grid (Triangulation<3> &coarse_grid) - { - GridGenerator::hyper_ball (coarse_grid); - static HyperBallBoundary<3> boundary; - coarse_grid.set_boundary (0, boundary); - coarse_grid.refine_global (1); - } - -} - -// @sect4{Discussion} -// -// As you have now read through this -// framework, you may be wondering -// why we have not chosen to -// implement the classes implementing -// a certain setup (like the -// ``CurvedRidges'' class) directly -// as classes derived from -// ``Data::SetUpBase''. Indeed, we -// could have done very well so. The -// only reason is that then we would -// have to have member variables for -// the solution and right hand side -// classes in the ``CurvedRidges'' -// class, as well as member functions -// overloading the abstract functions -// of the base class giving access to -// these member variables. The -// ``SetUp'' class has the sole -// reason to relieve us from the need -// to reiterate these member -// variables and functions that would -// be necessary in all such -// classes. In some way, the template -// mechanism here only provides a way -// to have default implementations -// for a number of functions that -// depend on external quantities and -// can thus not be provided using -// normal virtual functions, at least -// not without the help of templates. -// -// However, there might be good -// reasons to actually implement -// classes derived from -// ``Data::SetUpBase'', for example -// if the solution or right hand side -// classes require constructors that -// take arguments, which the -// ``Data::SetUpBase'' class cannot -// provide. In that case, subclassing -// is a worthwhile strategy. Other -// possibilities for special cases -// are to derive from -// ``Data::SetUp'' where -// ``SomeSetUp'' denotes a class, or -// even to explicitly specialize -// ``Data::SetUp''. The -// latter allows to transparently use -// the way the ``SetUp'' class is -// used for other set-ups, but with -// special actions taken for special -// arguments. -// -// A final observation favoring the -// approach taken here is the -// following: we have found numerous -// times that when starting a -// project, the number of parameters -// (usually boundary values, right -// hand side, coarse grid, just as -// here) was small, and the number of -// test cases was small as well. One -// then starts out by handcoding them -// into a number of ``switch'' -// statements. Over time, projects -// grow, and so does the number of -// test cases. The number of -// ``switch'' statements grows with -// that, and their length as well, -// and one starts to find ways to -// consider impossible examples where -// domains, boundary values, and -// right hand sides do not fit -// together any more, and starts -// loosing the overview over the -// whole structure. Encapsulating -// everything belonging to a certain -// test case into a structure of its -// own has proven worthwhile for -// this, as it keeps everything that -// belongs to one test case in one -// place. Furthermore, it allows to -// put these things all in one or -// more files that are only devoted -// to test cases and their data, -// without having to bring their -// actual implementation into contact -// with the rest of the program. - - -// @sect3{Dual functionals} - -// As with the other components of -// the program, we put everything we -// need to describe dual functionals -// into a namespace of its own, and -// define an abstract base class that -// provides the interface the class -// solving the dual problem needs for -// its work. -// -// We will then implement two such -// classes, for the evaluation of a -// point value and of the derivative -// of the solution at that point. For -// these functionals we already have -// the corresponding evaluation -// objects, so they are comlementary. -namespace DualFunctional -{ - // @sect4{The DualFunctionalBase class} - - // First start with the base class - // for dual functionals. Since for - // linear problems the - // characteristics of the dual - // problem play a role only in the - // right hand side, we only need to - // provide for a function that - // assembles the right hand side - // for a given discretization: - template - class DualFunctionalBase : public Subscriptor - { - public: - virtual - void - assemble_rhs (const DoFHandler &dof_handler, - Vector &rhs) const = 0; - }; - - - // @sect4{The PointValueEvaluation class} - - // As a first application, we - // consider the functional - // corresponding to the evaluation - // of the solution's value at a - // given point which again we - // assume to be a vertex. Apart - // from the constructor that takes - // and stores the evaluation point, - // this class consists only of the - // function that implements - // assembling the right hand side. - template - class PointValueEvaluation : public DualFunctionalBase - { - public: - PointValueEvaluation (const Point &evaluation_point); - - virtual - void - assemble_rhs (const DoFHandler &dof_handler, - Vector &rhs) const; - - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); - - - const Point evaluation_point; - }; - - - template - PointValueEvaluation:: - PointValueEvaluation (const Point &evaluation_point) - : - evaluation_point (evaluation_point) - {} - - - // As for doing the main purpose of - // the class, assembling the right - // hand side, let us first consider - // what is necessary: The right - // hand side of the dual problem is - // a vector of values J(phi_i), - // where J is the error functional, - // and phi_i is the i-th shape - // function. Here, J is the - // evaluation at the point x0, - // i.e. J(phi_i)=phi_i(x0). - // - // Now, we have assumed that the - // evaluation point is a - // vertex. Thus, for the usual - // finite elements we might be - // using in this program, we can - // take for granted that at such a - // point exactly one shape function - // is nonzero, and in particular - // has the value one. Thus, we set - // the right hand side vector to - // all-zeros, then seek for the - // shape function associated with - // that point and set the - // corresponding value of the right - // hand side vector to one: - template - void - PointValueEvaluation:: - assemble_rhs (const DoFHandler &dof_handler, - Vector &rhs) const - { - // So, first set everything to - // zeros... - rhs.reinit (dof_handler.n_dofs()); - - // ...then loop over cells and - // find the evaluation point - // among the vertices (or very - // close to a vertex, which may - // happen due to floating point - // round-off): - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - for (unsigned int vertex=0; - vertex::vertices_per_cell; - ++vertex) - if (cell->vertex(vertex).distance(evaluation_point) - < cell->diameter()*1e-8) - { - // Ok, found, so set - // corresponding entry, - // and leave function - // since we are finished: - rhs(cell->vertex_dof_index(vertex,0)) = 1; - return; - }; - - // Finally, a sanity check: if we - // somehow got here, then we must - // have missed the evaluation - // point, so raise an exception - // unconditionally: - AssertThrow (false, ExcEvaluationPointNotFound(evaluation_point)); - } - - - // @sect4{The PointXDerivativeEvaluation class} - - // As second application, we again - // consider the evaluation of the - // x-derivative of the solution at - // one point. Again, the - // declaration of the class, and - // the implementation of its - // constructor is not too - // interesting: - template - class PointXDerivativeEvaluation : public DualFunctionalBase - { - public: - PointXDerivativeEvaluation (const Point &evaluation_point); - - virtual - void - assemble_rhs (const DoFHandler &dof_handler, - Vector &rhs) const; - - DeclException1 (ExcEvaluationPointNotFound, - Point, - << "The evaluation point " << arg1 - << " was not found among the vertices of the present grid."); - - - const Point evaluation_point; - }; - - - template - PointXDerivativeEvaluation:: - PointXDerivativeEvaluation (const Point &evaluation_point) - : - evaluation_point (evaluation_point) - {} - - - // What is interesting is the - // implementation of this - // functional: here, - // J(phi_i)=d/dx phi_i(x0). - // - // We could, as in the - // implementation of the respective - // evaluation object take the - // average of the gradients of each - // shape function phi_i at this - // evaluation point. However, we - // take a slightly different - // approach: we simply take the - // average over all cells that - // surround this point. The - // question which cells - // ``surrounds'' the evaluation - // point is made dependent on the - // mesh width by including those - // cells for which the distance of - // the cell's midpoint to the - // evaluation point is less than - // the cell's diameter. - // - // Taking the average of the - // gradient over the area/volume of - // these cells leads to a dual - // solution which is very close to - // the one which would result from - // the point evaluation of the - // gradient. It is simple to - // justify theoretically that this - // does not change the method - // significantly. - template - void - PointXDerivativeEvaluation:: - assemble_rhs (const DoFHandler &dof_handler, - Vector &rhs) const - { - // Again, first set all entries - // to zero: - rhs.reinit (dof_handler.n_dofs()); - - // Initialize a ``FEValues'' - // object with a quadrature - // formula, have abbreviations - // for the number of quadrature - // points and shape functions... - QGauss quadrature(4); - MappingQ mapping (4); - FEValues fe_values (mapping, dof_handler.get_fe(), quadrature, - update_gradients | - update_q_points | - update_JxW_values); - const unsigned int n_q_points = fe_values.n_quadrature_points; - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; - - // ...and have two objects that - // are used to store the global - // indices of the degrees of - // freedom on a cell, and the - // values of the gradients of the - // shape functions at the - // quadrature points: - Vector cell_rhs (dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); - - // Finally have a variable in - // which we will sum up the - // area/volume of the cells over - // which we integrate, by - // integrating the unit functions - // on these cells: - double total_volume = 0; - - // Then start the loop over all - // cells, and select those cells - // which are close enough to the - // evaluation point: - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - if (cell->center().distance(evaluation_point) <= - cell->diameter()) - { - // If we have found such a - // cell, then initialize - // the ``FEValues'' object - // and integrate the - // x-component of the - // gradient of each shape - // function, as well as the - // unit function for the - // total area/volume. - fe_values.reinit (cell); - cell_rhs = 0; - - for (unsigned int q=0; qget_dof_indices (local_dof_indices); - for (unsigned int i=0; i 0, - ExcEvaluationPointNotFound(evaluation_point)); - - // Finally, we have by now only - // integrated the gradients of - // the shape functions, not - // taking their mean value. We - // fix this by dividing by the - // measure of the volume over - // which we have integrated: - rhs.scale (1./total_volume); - } - - -} - - -// @sect3{Extending the LaplaceSolver namespace} -namespace LaplaceSolver -{ - - // @sect4{The DualSolver class} - - // In the same way as the - // ``PrimalSolver'' class above, we - // now implement a - // ``DualSolver''. It has all the - // same features, the only - // difference is that it does not - // take a function object denoting - // a right hand side object, but - // now takes a - // ``DualFunctionalBase'' object - // that will assemble the right - // hand side vector of the dual - // problem. The rest of the class - // is rather trivial. - // - // Since both primal and dual - // solver will use the same - // triangulation, but different - // discretizations, it now becomes - // clear why we have made the - // ``Base'' class a virtual one: - // since the final class will be - // derived from both - // ``PrimalSolver'' as well as - // ``DualSolver'', it would have - // two ``Base'' instances, would we - // not have marked the inheritance - // as virtual. Since in many - // applications the base class - // would store much more - // information than just the - // triangulation which needs to be - // shared between primal and dual - // solvers, we do not usually want - // to use two such base classes. - template - class DualSolver : public Solver - { - public: - DualSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const DualFunctional::DualFunctionalBase &dual_functional); - - virtual - void - solve_problem (); - - virtual - unsigned int - n_dofs () const; - - virtual - void - postprocess (const Evaluation::EvaluationBase &postprocessor) const; - - - const SmartPointer > dual_functional; - virtual void assemble_rhs (Vector &rhs) const; - - static const ZeroFunction boundary_values; - - // Same as above -- make a - // derived class a friend of - // this one: - friend class WeightedResidual; - }; - - template - const ZeroFunction DualSolver::boundary_values; - - template - DualSolver:: - DualSolver (Triangulation &triangulation, - const FiniteElement &fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const DualFunctional::DualFunctionalBase &dual_functional) - : - Base (triangulation), - Solver (triangulation, fe, - quadrature, face_quadrature, - boundary_values), - dual_functional (&dual_functional) - {} - - - template - void - DualSolver::solve_problem () - { - Solver::solve_problem (); - } - - - - template - unsigned int - DualSolver::n_dofs() const - { - return Solver::n_dofs(); - } - - - template - void - DualSolver:: - postprocess (const Evaluation::EvaluationBase &postprocessor) const - { - Solver::postprocess(postprocessor); - } - - - - template - void - DualSolver:: - assemble_rhs (Vector &rhs) const - { - dual_functional->assemble_rhs (this->dof_handler, rhs); - } - - - // @sect4{The WeightedResidual class} - - // Here finally comes the main - // class of this program, the one - // that implements the dual - // weighted residual error - // estimator. It joins the primal - // and dual solver classes to use - // them for the computation of - // primal and dual solutions, and - // implements the error - // representation formula for use - // as error estimate and mesh - // refinement. - // - // The first few of the functions - // of this class are mostly - // overriders of the respective - // functions of the base class: - template - class WeightedResidual : public PrimalSolver, - public DualSolver - { - public: - WeightedResidual (Triangulation &coarse_grid, - const FiniteElement &primal_fe, - const FiniteElement &dual_fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &boundary_values, - const DualFunctional::DualFunctionalBase &dual_functional); - - virtual - void - solve_problem (); - - virtual - void - postprocess (const Evaluation::EvaluationBase &postprocessor) const; - - virtual - unsigned int - n_dofs () const; - - virtual void refine_grid (); - - virtual - void - output_solution () const; - - - // In the private section, we - // have two functions that are - // used to call the - // ``solve_problem'' functions - // of the primal and dual base - // classes. These two functions - // will be called in parallel - // by the ``solve_problem'' - // function of this class. - void solve_primal_problem (); - void solve_dual_problem (); - // Then declare abbreviations - // for active cell iterators, - // to avoid that we have to - // write this lengthy name - // over and over again: - - typedef - typename DoFHandler::active_cell_iterator - active_cell_iterator; - - // Next, declare a data type - // that we will us to store the - // contribution of faces to the - // error estimator. The idea is - // that we can compute the face - // terms from each of the two - // cells to this face, as they - // are the same when viewed - // from both sides. What we - // will do is to compute them - // only once, based on some - // rules explained below which - // of the two adjacent cells - // will be in charge to do - // so. We then store the - // contribution of each face in - // a map mapping faces to their - // values, and only collect the - // contributions for each cell - // by looping over the cells a - // second time and grabbing the - // values from the map. - // - // The data type of this map is - // declared here: - typedef - typename std::map::face_iterator,double> - FaceIntegrals; - - // In the computation of the - // error estimates on cells and - // faces, we need a number of - // helper objects, such as - // ``FEValues'' and - // ``FEFaceValues'' functions, - // but also temporary objects - // storing the values and - // gradients of primal and dual - // solutions, for - // example. These fields are - // needed in the three - // functions that do the - // integration on cells, and - // regular and irregular faces, - // respectively. - // - // There are three reasonable - // ways to provide these - // fields: first, as local - // variables in the function - // that needs them; second, as - // member variables of this - // class; third, as arguments - // passed to that function. - // - // These three alternatives all - // have drawbacks: the third - // that their number is not - // neglectable and would make - // calling these functions a - // lengthy enterprise. The - // second has the drawback that - // it disallows - // parallelization, since the - // threads that will compute - // the error estimate have to - // have their own copies of - // these variables each, so - // member variables of the - // enclosing class will not - // work. The first approach, - // although straightforward, - // has a subtle but important - // drawback: we will call these - // functions over and over - // again, many thousand times - // maybe; it has now turned out - // that allocating vectors and - // other objects that need - // memory from the heap is an - // expensive business in terms - // of run-time, since memory - // allocation is expensive when - // several threads are - // involved. In our experience, - // more than 20 per cent of the - // total run time of error - // estimation functions are due - // to memory allocation, if - // done on a per-call level. It - // is thus significantly better - // to allocate the memory only - // once, and recycle the - // objects as often as - // possible. - // - // What to do? Our answer is to - // use a variant of the third - // strategy, namely generating - // these variables once in the - // main function of each - // thread, and passing them - // down to the functions that - // do the actual work. To avoid - // that we have to give these - // functions a dozen or so - // arguments, we pack all these - // variables into two - // structures, one which is - // used for the computations on - // cells, the other doing them - // on the faces. Instead of - // many individual objects, we - // will then only pass one such - // object to these functions, - // making their calling - // sequence simpler. - struct CellData - { - MappingQ mapping; - FEValues fe_values; - const SmartPointer > right_hand_side; - - std::vector cell_residual; - std::vector rhs_values; - std::vector dual_weights; - typename std::vector > cell_grad_grads; - CellData (const FiniteElement &fe, - const Quadrature &quadrature, - const Function &right_hand_side); - }; - - struct FaceData - { - MappingQ mapping; - FEFaceValues fe_face_values_cell; - FEFaceValues fe_face_values_neighbor; - FESubfaceValues fe_subface_values_cell; - - std::vector jump_residual; - std::vector dual_weights; - typename std::vector > cell_grads; - typename std::vector > neighbor_grads; - FaceData (const FiniteElement &fe, - const Quadrature &face_quadrature); - }; - - - - // Regarding the evaluation of - // the error estimator, we have - // two driver functions that do - // this: the first is called to - // generate the cell-wise - // estimates, and splits up the - // task in a number of threads - // each of which work on a - // subset of the cells. The - // first function will run the - // second for each of these - // threads: - void estimate_error (Vector &error_indicators) const; - - void estimate_some (const Vector &primal_solution, - const Vector &dual_weights, - const unsigned int n_threads, - const unsigned int this_thread, - Vector &error_indicators, - FaceIntegrals &face_integrals) const; - - // Then we have functions that - // do the actual integration of - // the error representation - // formula. They will treat the - // terms on the cell interiors, - // on those faces that have no - // hanging nodes, and on those - // faces with hanging nodes, - // respectively: - void - integrate_over_cell (const active_cell_iterator &cell, - const unsigned int cell_index, - const Vector &primal_solution, - const Vector &dual_weights, - CellData &cell_data, - Vector &error_indicators) const; - - void - integrate_over_regular_face (const active_cell_iterator &cell, - const unsigned int face_no, - const Vector &primal_solution, - const Vector &dual_weights, - FaceData &face_data, - FaceIntegrals &face_integrals) const; - void - integrate_over_irregular_face (const active_cell_iterator &cell, - const unsigned int face_no, - const Vector &primal_solution, - const Vector &dual_weights, - FaceData &face_data, - FaceIntegrals &face_integrals) const; - }; - - - - // In the implementation of this - // class, we first have the - // constructors of the ``CellData'' - // and ``FaceData'' member classes, - // and the ``WeightedResidual'' - // constructor. They only - // initialize fields to their - // correct lengths, so we do not - // have to discuss them to length. - template - WeightedResidual::CellData:: - CellData (const FiniteElement &fe, - const Quadrature &quadrature, - const Function &right_hand_side) - : - mapping (4), - fe_values (mapping, fe, quadrature, - update_values | - update_second_derivatives | - update_q_points | - update_JxW_values), - right_hand_side (&right_hand_side) - { - const unsigned int n_q_points - = quadrature.n_quadrature_points; - - cell_residual.resize(n_q_points); - rhs_values.resize(n_q_points); - dual_weights.resize(n_q_points); - cell_grad_grads.resize(n_q_points); - } - - - - template - WeightedResidual::FaceData:: - FaceData (const FiniteElement &fe, - const Quadrature &face_quadrature) - : - mapping (4), - fe_face_values_cell (mapping, fe, face_quadrature, - update_values | - update_gradients | - update_JxW_values | - update_normal_vectors), - fe_face_values_neighbor (mapping, fe, face_quadrature, - update_values | - update_gradients | - update_JxW_values | - update_normal_vectors), - fe_subface_values_cell (mapping, fe, face_quadrature, - update_gradients) - { - const unsigned int n_face_q_points - = face_quadrature.n_quadrature_points; - - jump_residual.resize(n_face_q_points); - dual_weights.resize(n_face_q_points); - cell_grads.resize(n_face_q_points); - neighbor_grads.resize(n_face_q_points); - } - - - - - template - WeightedResidual:: - WeightedResidual (Triangulation &coarse_grid, - const FiniteElement &primal_fe, - const FiniteElement &dual_fe, - const Quadrature &quadrature, - const Quadrature &face_quadrature, - const Function &rhs_function, - const Function &bv, - const DualFunctional::DualFunctionalBase &dual_functional) - : - Base (coarse_grid), - PrimalSolver (coarse_grid, primal_fe, - quadrature, face_quadrature, - rhs_function, bv), - DualSolver (coarse_grid, dual_fe, - quadrature, face_quadrature, - dual_functional) - {} - - - // The next five functions are - // boring, as they simply relay - // their work to the base - // classes. The first calls the - // primal and dual solvers in - // parallel, while postprocessing - // the solution and retrieving the - // number of degrees of freedom is - // done by the primal class. - template - void - WeightedResidual::solve_problem () - { - Threads::ThreadGroup<> threads; - threads += Threads::new_thread (&WeightedResidual::solve_primal_problem, *this); - threads += Threads::new_thread (&WeightedResidual::solve_dual_problem,*this); - threads.join_all (); - } - - - template - void - WeightedResidual::solve_primal_problem () - { - PrimalSolver::solve_problem (); - } - - template - void - WeightedResidual::solve_dual_problem () - { - DualSolver::solve_problem (); - } - - - template - void - WeightedResidual:: - postprocess (const Evaluation::EvaluationBase &postprocessor) const - { - PrimalSolver::postprocess (postprocessor); - } - - - template - unsigned int - WeightedResidual::n_dofs () const - { - return PrimalSolver::n_dofs(); - } - - - - // Now, it is becoming more - // interesting: the ``refine_grid'' - // function asks the error - // estimator to compute the - // cell-wise error indicators, then - // uses their absolute values for - // mesh refinement. - template - void - WeightedResidual::refine_grid () - { - // First call the function that - // computes the cell-wise and - // global error: - Vector error_indicators (this->triangulation->n_active_cells()); - estimate_error (error_indicators); - - // Then note that marking cells - // for refinement or coarsening - // only works if all indicators - // are positive, to allow their - // comparison. Thus, drop the - // signs on all these indicators: - for (Vector::iterator i=error_indicators.begin(); - i != error_indicators.end(); ++i) - *i = std::fabs (*i); - - // Finally, we can select between - // different strategies for - // refinement. The default here - // is to refine those cells with - // the largest error indicators - // that make up for a total of 80 - // per cent of the error, while - // we coarsen those with the - // smallest indicators that make - // up for the bottom 2 per cent - // of the error. - GridRefinement::refine_and_coarsen_fixed_fraction (*this->triangulation, - error_indicators, - 0.5/std::sqrt(1.+step), - 0.2/std::sqrt(1.+step)); - this->triangulation->execute_coarsening_and_refinement (); - } - - - // Since we want to output both the - // primal and the dual solution, we - // overload the ``output_solution'' - // function. The only interesting - // feature of this function is that - // the primal and dual solutions - // are defined on different finite - // element spaces, which is not the - // format the ``DataOut'' class - // expects. Thus, we have to - // transfer them to a common finite - // element space. Since we want the - // solutions only to see them - // qualitatively, we contend - // ourselves with interpolating the - // dual solution to the (smaller) - // primal space. For the - // interpolation, there is a - // library function, that takes a - // ``ConstraintMatrix'' object - // including the hanging node - // constraints. The rest is - // standard. - // - // There is, however, one - // work-around worth mentioning: in - // this function, as in a couple of - // following ones, we have to - // access the ``DoFHandler'' - // objects and solutions of both - // the primal as well as of the - // dual solver. Since these are - // members of the ``Solver'' base - // class which exists twice in the - // class hierarchy leading to the - // present class (once as base - // class of the ``PrimalSolver'' - // class, once as base class of the - // ``DualSolver'' class), we have - // to disambiguate accesses to them - // by telling the compiler a member - // of which of these two instances - // we want to access. The way to do - // this would be identify the - // member by pointing a path - // through the class hierarchy - // which disambiguates the base - // class, for example writing - // ``PrimalSolver::dof_handler'' to - // denote the member variable - // ``dof_handler'' from the - // ``Solver'' base class of the - // ``PrimalSolver'' - // class. Unfortunately, this - // confuses gcc's version 2.96 (a - // version that was intended as a - // development snapshot, but - // delivered as system compiler by - // Red Hat in their 7.x releases) - // so much that it bails out and - // refuses to compile the code. - // - // Thus, we have to work around - // this problem. We do this by - // introducing references to the - // ``PrimalSolver'' and - // ``DualSolver'' components of the - // ``WeightedResidual'' object at - // the beginning of the - // function. Since each of these - // has an unambiguous base class - // ``Solver'', we can access the - // member variables we want through - // these references. However, we - // are now accessing protected - // member variables of these - // classes through a pointer other - // than the ``this'' pointer (in - // fact, this is of course the - // ``this'' pointer, but not - // explicitly). This finally is the - // reason why we had to declare the - // present class a friend of the - // classes we so access. - template - void - WeightedResidual::output_solution () const - { - const PrimalSolver &primal_solver = *this; - const DualSolver &dual_solver = *this; - - ConstraintMatrix primal_hanging_node_constraints; - DoFTools::make_hanging_node_constraints (primal_solver.dof_handler, - primal_hanging_node_constraints); - primal_hanging_node_constraints.close(); - Vector dual_solution (primal_solver.dof_handler.n_dofs()); - FETools::interpolate (dual_solver.dof_handler, - dual_solver.solution, - primal_solver.dof_handler, - primal_hanging_node_constraints, - dual_solution); - - // approximate error, gradient, - // and second derivative - // information as cell information - Vector error_indicators (this->triangulation->n_active_cells()); - Vector gradient_indicators (this->triangulation->n_active_cells()); - Vector second_indicators (this->triangulation->n_active_cells()); - { - MappingQ mapping(4); - KellyErrorEstimator::estimate (mapping, primal_solver.dof_handler, - QGauss(3), - typename FunctionMap::type(), - primal_solver.solution, - error_indicators); - - DerivativeApproximation:: - approximate_gradient (mapping, - primal_solver.dof_handler, - primal_solver.solution, - gradient_indicators); - - DerivativeApproximation:: - approximate_second_derivative (mapping, - primal_solver.dof_handler, - primal_solver.solution, - second_indicators); - - } - // distribute cell to dof vectors - Vector x_error_indicators (primal_solver.dof_handler.n_dofs()); - Vector x_gradient_indicators (primal_solver.dof_handler.n_dofs()); - Vector x_second_indicators (primal_solver.dof_handler.n_dofs()); - DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler, - error_indicators, - x_error_indicators); - DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler, - gradient_indicators, - x_gradient_indicators); - DoFTools::distribute_cell_to_dof_vector (primal_solver.dof_handler, - second_indicators, - x_second_indicators); - - - - // we generate too much output in - // 3d. instead of doing it that - // way, simply generate a coarser - // mesh and output from there - Triangulation coarser_mesh; - coarser_mesh.copy_triangulation (*this->triangulation); - for (typename Triangulation::active_cell_iterator - cell = coarser_mesh.begin_active(); - cell != coarser_mesh.end(); ++cell) - cell->set_coarsen_flag(); - coarser_mesh.execute_coarsening_and_refinement (); - - // next generate a DoF handler on - // that mesh and a map fron one - // to the other mesh - DoFHandler coarser_dof_handler (coarser_mesh); - coarser_dof_handler.distribute_dofs (primal_solver.dof_handler.get_fe()); - InterGridMap > coarse_to_fine_map; - coarse_to_fine_map.make_mapping (coarser_dof_handler, - primal_solver.dof_handler); - - // finally we have to transfer - // the data vectors - Vector coarse_primal_solution (coarser_dof_handler.n_dofs()); - Vector coarse_dual_solution (coarser_dof_handler.n_dofs()); - Vector coarse_error_indicators (coarser_dof_handler.n_dofs()); - Vector coarse_gradient_indicators (coarser_dof_handler.n_dofs()); - Vector coarse_second_indicators (coarser_dof_handler.n_dofs()); - - Vector tmp (coarser_dof_handler.get_fe().dofs_per_cell); - for (typename DoFHandler::active_cell_iterator - cell = coarser_dof_handler.begin_active(); - cell != coarser_dof_handler.end(); ++cell) - { - coarse_to_fine_map[cell]->get_interpolated_dof_values (primal_solver.solution,tmp); - cell->set_dof_values (tmp, coarse_primal_solution); - - coarse_to_fine_map[cell]->get_interpolated_dof_values (dual_solution,tmp); - cell->set_dof_values (tmp, coarse_dual_solution); - - coarse_to_fine_map[cell]->get_interpolated_dof_values (x_error_indicators,tmp); - cell->set_dof_values (tmp, coarse_error_indicators); - - coarse_to_fine_map[cell]->get_interpolated_dof_values (x_gradient_indicators,tmp); - cell->set_dof_values (tmp, coarse_gradient_indicators); - - coarse_to_fine_map[cell]->get_interpolated_dof_values (x_second_indicators,tmp); - cell->set_dof_values (tmp, coarse_second_indicators); - } - - { - DataOut data_out; - data_out.attach_dof_handler (coarser_dof_handler); - data_out.add_data_vector (coarse_primal_solution, "primal_solution"); - data_out.add_data_vector (coarse_dual_solution, "dual_solution"); - data_out.add_data_vector (coarse_error_indicators, "errors"); - data_out.add_data_vector (coarse_gradient_indicators, "gradient"); - data_out.add_data_vector (coarse_second_indicators, "second_derivatives"); - data_out.build_patches (); - -#ifdef HAVE_STD_STRINGSTREAM - std::ostringstream filename; -#else - std::ostrstream filename; -#endif - filename << "spec2006-447.dealII/" - << "solution-" - << this->refinement_cycle - << ".gmv" - << std::ends; -#ifdef HAVE_STD_STRINGSTREAM - std::ofstream out (filename.str().c_str()); -#else - std::ofstream out (filename.str()); -#endif - - data_out.write_gmv (out); - } - - } - - - // @sect3{Estimating errors} - - // @sect4{Error estimation driver functions} - // - // As for the actual computation of - // error estimates, let's start - // with the function that drives - // all this, i.e. calls those - // functions that actually do the - // work, and finally collects the - // results. - - template - void - WeightedResidual:: - estimate_error (Vector &error_indicators) const - { - const PrimalSolver &primal_solver = *this; - const DualSolver &dual_solver = *this; - - // The first task in computing - // the error is to set up vectors - // that denote the primal - // solution, and the weights - // (z-z_h)=(z-I_hz), both in the - // finite element space for which - // we have computed the dual - // solution. For this, we have to - // interpolate the primal - // solution to the dual finite - // element space, and to subtract - // the interpolation of the - // computed dual solution to the - // primal finite element - // space. Fortunately, the - // library provides functions for - // the interpolation into larger - // or smaller finite element - // spaces, so this is mostly - // obvious. - // - // First, let's do that for the - // primal solution: it is - // cell-wise interpolated into - // the finite element space in - // which we have solved the dual - // problem: But, again as in the - // ``WeightedResidual::output_solution'' - // function we first need to - // create a ConstraintMatrix - // including the hanging node - // constraints, but this time of - // the dual finite element space. - ConstraintMatrix dual_hanging_node_constraints; - DoFTools::make_hanging_node_constraints (dual_solver.dof_handler, - dual_hanging_node_constraints); - dual_hanging_node_constraints.close(); - Vector primal_solution (dual_solver.dof_handler.n_dofs()); - FETools::interpolate (primal_solver.dof_handler, - primal_solver.solution, - dual_solver.dof_handler, - dual_hanging_node_constraints, - primal_solution); - - // Then for computing the - // interpolation of the - // numerically approximated dual - // solution z into the finite - // element space of the primal - // solution and subtracting it - // from z: use the - // ``interpolate_difference'' - // function, that gives (z-I_hz) - // in the element space of the - // dual solution. - ConstraintMatrix primal_hanging_node_constraints; - DoFTools::make_hanging_node_constraints (primal_solver.dof_handler, - primal_hanging_node_constraints); - primal_hanging_node_constraints.close(); - Vector dual_weights (dual_solver.dof_handler.n_dofs()); - FETools::interpolation_difference (dual_solver.dof_handler, - dual_hanging_node_constraints, - dual_solver.solution, - primal_solver.dof_handler, - primal_hanging_node_constraints, - dual_weights); - - // Note that this could probably - // have been more efficient since - // those constraints have been - // used previously when - // assembling matrix and right - // hand side for the primal - // problem and writing out the - // dual solution. We leave the - // optimization of the program in - // this respect as an exercise. - - // Having computed the dual - // weights we now proceed with - // computing the cell and face - // residuals of the primal - // solution. First we set up a - // map between face iterators and - // their jump term contributions - // of faces to the error - // estimator. The reason is that - // we compute the jump terms only - // once, from one side of the - // face, and want to collect them - // only afterwards when looping - // over all cells a second time. - // - // We initialize this map already - // with a value of -1e20 for all - // faces, since this value will - // strike in the results if - // something should go wrong and - // we fail to compute the value - // for a face for some - // reason. Secondly, we - // initialize the map once before - // we branch to different threads - // since this way the map's - // structure is no more modified - // by the individual threads, - // only existing entries are set - // to new values. This relieves - // us from the necessity to - // synchronise the threads - // through a mutex each time they - // write to (and modify the - // structure of) this map. - FaceIntegrals face_integrals; - for (active_cell_iterator cell=dual_solver.dof_handler.begin_active(); - cell!=dual_solver.dof_handler.end(); - ++cell) - for (unsigned int face_no=0; - face_no::faces_per_cell; - ++face_no) - face_integrals[cell->face(face_no)] = -1e20; - - // Then set up a vector with - // error indicators. Reserve one - // slot for each cell and set it - // to zero. - error_indicators.reinit (dual_solver.dof_handler - .get_tria().n_active_cells()); - - // Now start a number of threads - // which compute the error - // formula on parts of all the - // cells, and once they are all - // started wait until they have - // all finished: - const unsigned int n_threads = multithread_info.n_default_threads; - Threads::ThreadGroup<> threads; - for (unsigned int i=0; i::estimate_some, *this, - primal_solution, - dual_weights, - n_threads, i, - error_indicators, - face_integrals); - threads.join_all(); - - // Once the error contributions - // are computed, sum them up. For - // this, note that the cell terms - // are already set, and that only - // the edge terms need to be - // collected. Thus, loop over all - // cells and their faces, make - // sure that the contributions of - // each of the faces are there, - // and add them up. Only take - // minus one half of the jump - // term, since the other half - // will be taken by the - // neighboring cell. - unsigned int present_cell=0; - for (active_cell_iterator cell=dual_solver.dof_handler.begin_active(); - cell!=dual_solver.dof_handler.end(); - ++cell, ++present_cell) - for (unsigned int face_no=0; face_no::faces_per_cell; - ++face_no) - { - Assert(face_integrals.find(cell->face(face_no)) != - face_integrals.end(), - ExcInternalError()); - error_indicators(present_cell) - -= 0.5*face_integrals[cell->face(face_no)]; - }; - std::cout << " Estimated error=" - << std::accumulate (error_indicators.begin(), - error_indicators.end(), 0.) - << std::endl; - } - - - // @sect4{Estimating on a subset of cells} - - // Next we have the function that - // is called to estimate the error - // on a subset of cells. The - // function may be called multiply - // if the library was configured to - // use multi-threading. Here it - // goes: - template - void - WeightedResidual:: - estimate_some (const Vector &primal_solution, - const Vector &dual_weights, - const unsigned int n_threads, - const unsigned int this_thread, - Vector &error_indicators, - FaceIntegrals &face_integrals) const - { - const PrimalSolver &primal_solver = *this; - const DualSolver &dual_solver = *this; - - // At the beginning, we - // initialize two variables for - // each thread which may be - // running this function. The - // reason for these functions was - // discussed above, when the - // respective classes were - // discussed, so we here only - // point out that since they are - // local to the function that is - // spawned when running more than - // one thread, the data of these - // objects exists actually once - // per thread, so we don't have - // to take care about - // synchronising access to them. - CellData cell_data (*dual_solver.fe, - *dual_solver.quadrature, - *primal_solver.rhs_function); - FaceData face_data (*dual_solver.fe, - *dual_solver.face_quadrature); - - // Then calculate the start cell - // for this thread. We let the - // different threads run on - // interleaved cells, i.e. for - // example if we have 4 threads, - // then the first thread treates - // cells 0, 4, 8, etc, while the - // second threads works on cells 1, - // 5, 9, and so on. The reason is - // that it takes vastly more time - // to work on cells with hanging - // nodes than on regular cells, but - // such cells are not evenly - // distributed across the range of - // cell iterators, so in order to - // have the different threads do - // approximately the same amount of - // work, we have to let them work - // interleaved to the effect of a - // pseudorandom distribution of the - // `hard' cells to the different - // threads. - active_cell_iterator cell=dual_solver.dof_handler.begin_active(); - for (unsigned int t=0; - (t::faces_per_cell; - ++face_no) - { - // First, if this face is - // part of the boundary, - // then there is nothing - // to do. However, to - // make things easier - // when summing up the - // contributions of the - // faces of cells, we - // enter this face into - // the list of faces with - // a zero contribution to - // the error. - if (cell->face(face_no)->at_boundary()) - { - face_integrals[cell->face(face_no)] = 0; - continue; - }; - - // Next, note that since - // we want to compute the - // jump terms on each - // face only once - // although we access it - // twice (if it is not at - // the boundary), we have - // to define some rules - // who is responsible for - // computing on a face: - // - // First, if the - // neighboring cell is on - // the same level as this - // one, i.e. neither - // further refined not - // coarser, then the one - // with the lower index - // within this level does - // the work. In other - // words: if the other - // one has a lower index, - // then skip work on this - // face: - if ((cell->neighbor(face_no)->has_children() == false) && - (cell->neighbor(face_no)->level() == cell->level()) && - (cell->neighbor(face_no)->index() < cell->index())) - continue; - - // Likewise, we always - // work from the coarser - // cell if this and its - // neighbor differ in - // refinement. Thus, if - // the neighboring cell - // is less refined than - // the present one, then - // do nothing since we - // integrate over the - // subfaces when we visit - // the coarse cell. - if (cell->at_boundary(face_no) == false) - if (cell->neighbor(face_no)->level() < cell->level()) - continue; - - - // Now we know that we - // are in charge here, so - // actually compute the - // face jump terms. If - // the face is a regular - // one, i.e. the other - // side's cell is neither - // coarser not finer than - // this cell, then call - // one function, and if - // the cell on the other - // side is further - // refined, then use - // another function. Note - // that the case that the - // cell on the other side - // is coarser cannot - // happen since we have - // decided above that we - // handle this case when - // we pass over that - // other cell. - if (cell->face(face_no)->has_children() == false) - integrate_over_regular_face (cell, face_no, - primal_solution, - dual_weights, - face_data, - face_integrals); - else - integrate_over_irregular_face (cell, face_no, - primal_solution, - dual_weights, - face_data, - face_integrals); - }; - - // After computing the cell - // contributions and looping - // over the faces, go to the - // next cell for this - // thread. Note again that - // the cells for each of the - // threads are interleaved. - // If we are at the end of - // our workload, jump out - // of the loop. - for (unsigned int t=0; - ((t - void WeightedResidual:: - integrate_over_cell (const active_cell_iterator &cell, - const unsigned int cell_index, - const Vector &primal_solution, - const Vector &dual_weights, - CellData &cell_data, - Vector &error_indicators) const - { - // The tasks to be done are what - // appears natural from looking - // at the error estimation - // formula: first compute the the - // right hand side and the - // Laplacian of the numerical - // solution at the quadrature - // points for the cell residual, - cell_data.fe_values.reinit (cell); - cell_data.right_hand_side - ->value_list (cell_data.fe_values.get_quadrature_points(), - cell_data.rhs_values); - cell_data.fe_values.get_function_2nd_derivatives (primal_solution, - cell_data.cell_grad_grads); - - // ...then get the dual weights... - cell_data.fe_values.get_function_values (dual_weights, - cell_data.dual_weights); - - // ...and finally build the sum - // over all quadrature points and - // store it with the present - // cell: - double sum = 0; - for (unsigned int p=0; p - void WeightedResidual:: - integrate_over_regular_face (const active_cell_iterator &cell, - const unsigned int face_no, - const Vector &primal_solution, - const Vector &dual_weights, - FaceData &face_data, - FaceIntegrals &face_integrals) const - { - const unsigned int - n_q_points = face_data.fe_face_values_cell.n_quadrature_points; - - // The first step is to get the - // values of the gradients at the - // quadrature points of the - // finite element field on the - // present cell. For this, - // initialize the - // ``FEFaceValues'' object - // corresponding to this side of - // the face, and extract the - // gradients using that - // object. - face_data.fe_face_values_cell.reinit (cell, face_no); - face_data.fe_face_values_cell.get_function_grads (primal_solution, - face_data.cell_grads); - - // The second step is then to - // extract the gradients of the - // finite element solution at the - // quadrature points on the other - // side of the face, i.e. from - // the neighboring cell. - // - // For this, do a sanity check - // before: make sure that the - // neigbor actually exists (yes, - // we should not have come here - // if the neighbor did not exist, - // but in complicated software - // there are bugs, so better - // check this), and if this is - // not the case throw an error. - Assert (cell->neighbor(face_no).state() == IteratorState::valid, - ExcInternalError()); - // If we have that, then we need - // to find out with which face of - // the neighboring cell we have - // to work, i.e. the - // ``home-many''the neighbor the - // present cell is of the cell - // behind the present face. For - // this, there is a function, and - // we put the result into a - // variable with the name - // ``neighbor_neighbor'': - const unsigned int - neighbor_neighbor = cell->neighbor_of_neighbor (face_no); - // Then define an abbreviation - // for the neigbor cell, - // initialize the - // ``FEFaceValues'' object on - // that cell, and extract the - // gradients on that cell: - const active_cell_iterator neighbor = cell->neighbor(face_no); - face_data.fe_face_values_neighbor.reinit (neighbor, neighbor_neighbor); - face_data.fe_face_values_neighbor.get_function_grads (primal_solution, - face_data.neighbor_grads); - - // Now that we have the gradients - // on this and the neighboring - // cell, compute the jump - // residual by multiplying the - // jump in the gradient with the - // normal vector: - for (unsigned int p=0; pface(face_no)) != face_integrals.end(), - ExcInternalError()); - Assert (face_integrals[cell->face(face_no)] == -1e20, - ExcInternalError()); - - // ...then store computed value - // at assigned location. Note - // that the stored value does not - // contain the factor 1/2 that - // appears in the error - // representation. The reason is - // that the term actually does - // not have this factor if we - // loop over all faces in the - // triangulation, but only - // appears if we write it as a - // sum over all cells and all - // faces of each cell; we thus - // visit the same face twice. We - // take account of this by using - // this factor -1/2 later, when we - // sum up the contributions for - // each cell individually. - face_integrals[cell->face(face_no)] = face_integral; - } - - - // @sect4{Computing edge term error contributions - 2} - - // We are still missing the case of - // faces with hanging nodes. This - // is what is covered in this - // function: - template - void WeightedResidual:: - integrate_over_irregular_face (const active_cell_iterator &cell, - const unsigned int face_no, - const Vector &primal_solution, - const Vector &dual_weights, - FaceData &face_data, - FaceIntegrals &face_integrals) const - { - // First again two abbreviations, - // and some consistency checks - // whether the function is called - // only on faces for which it is - // supposed to be called: - const unsigned int - n_q_points = face_data.fe_face_values_cell.n_quadrature_points; - - const typename DoFHandler::cell_iterator - neighbor = cell->neighbor(face_no); - Assert (neighbor.state() == IteratorState::valid, - ExcInternalError()); - Assert (neighbor->has_children(), - ExcInternalError()); - - // Then find out which neighbor - // the present cell is of the - // adjacent cell. Note that we - // will operator on the children - // of this adjacent cell, but - // that their orientation is the - // same as that of their mother, - // i.e. the neigbor direction is - // the same. - const unsigned int - neighbor_neighbor = cell->neighbor_of_neighbor (face_no); - - // Then simply do everything we - // did in the previous function - // for one face for all the - // sub-faces now: - for (unsigned int subface_no=0; - subface_no::subfaces_per_face; - ++subface_no) - { - const active_cell_iterator neighbor_child - = cell->neighbor_child_on_subface (face_no, subface_no); - - // Now start the work by - // again getting the gradient - // of the solution first at - // this side of the - // interface, - face_data.fe_subface_values_cell.reinit (cell, face_no, subface_no); - face_data.fe_subface_values_cell.get_function_grads (primal_solution, - face_data.cell_grads); - // then at the other side, - face_data.fe_face_values_neighbor.reinit (neighbor_child, - neighbor_neighbor); - face_data.fe_face_values_neighbor.get_function_grads (primal_solution, - face_data.neighbor_grads); - - // and finally building the - // jump residuals. Since we - // take the normal vector - // from the other cell this - // time, revert the sign of - // the first term compared to - // the other function: - for (unsigned int p=0; pface(neighbor_neighbor)] - = face_integral; - }; - - // Once the contributions of all - // sub-faces are computed, loop - // over all sub-faces to collect - // and store them with the mother - // face for simple use when later - // collecting the error terms of - // cells. Again make safety - // checks that the entries for - // the sub-faces have been - // computed and do not carry an - // invalid value. - double sum = 0; - typename DoFHandler::face_iterator face = cell->face(face_no); - for (unsigned int subface_no=0; - subface_no::subfaces_per_face; - ++subface_no) - { - Assert (face_integrals.find(face->child(subface_no)) != - face_integrals.end(), - ExcInternalError()); - Assert (face_integrals[face->child(subface_no)] != -1e20, - ExcInternalError()); - - sum += face_integrals[face->child(subface_no)]; - }; - // Finally store the value with - // the parent face. - face_integrals[face] = sum; - } - -} - - -// @sect3{A simulation framework} - -// In the previous example program, -// we have had two functions that -// were used to drive the process of -// solving on subsequently finer -// grids. We extend this here to -// allow for a number of parameters -// to be passed to these functions, -// and put all of that into framework -// class. -// -// You will have noted that this -// program is built up of a number of -// small parts (evaluation functions, -// solver classes implementing -// various refinement methods, -// different dual functionals, -// different problem and data -// descriptions), which makes the -// program relatively simple to -// extend, but also allows to solve a -// large number of different problems -// by replacing one part by -// another. We reflect this -// flexibility by declaring a -// structure in the following -// framework class that holds a -// number of parameters that may be -// set to test various combinations -// of the parts of this program, and -// which can be used to test it at -// various problems and -// discretizations in a simple way. -template -struct Framework -{ -public: - // First, we declare two - // abbreviations for simple use - // of the respective data types: - typedef Evaluation::EvaluationBase Evaluator; - typedef std::list EvaluatorList; - - - // Then we have the structure - // which declares all the - // parameters that may be set. In - // the default constructor of the - // structure, these values are - // all set to default values, for - // simple use. - struct ProblemDescription - { - // First allow for the - // degrees of the piecewise - // polynomials by which the - // primal and dual problems - // will be discretized. They - // default to (bi-, - // tri-)linear ansatz - // functions for the primal, - // and (bi-, tri-)quadratic - // ones for the dual - // problem. If a refinement - // criterion is chosen that - // does not need the solution - // of a dual problem, the - // value of the dual finite - // element degree is of - // course ignored. - unsigned int primal_fe_degree; - unsigned int dual_fe_degree; - - // Then have an object that - // describes the problem - // type, i.e. right hand - // side, domain, boundary - // values, etc. The pointer - // needed here defaults to - // the Null pointer, i.e. you - // will have to set it in - // actual instances of this - // object to make it useful. - SmartPointer > data; - - // Since we allow to use - // different refinement - // criteria (global - // refinement, refinement by - // the Kelly error indicator, - // possibly with a weight, - // and using the dual - // estimator), define a - // number of enumeration - // values, and subsequently a - // variable of that type. It - // will default to - // ``dual_weighted_error_estimator''. - enum RefinementCriterion - { - dual_weighted_error_estimator, - global_refinement, - kelly_indicator, - weighted_kelly_indicator - }; - - RefinementCriterion refinement_criterion; - - // Next, an object that - // describes the dual - // functional. It is only - // needed if the dual - // weighted residual - // refinement is chosen, and - // also defaults to a Null - // pointer. - SmartPointer > dual_functional; - - // Then a list of evaluation - // objects. Its default value - // is empty, i.e. no - // evaluation objects. - EvaluatorList evaluator_list; - - // Next to last, a function - // that is used as a weight - // to the - // ``RefinementWeightedKelly'' - // class. The default value - // of this pointer is zero, - // but you have to set it to - // some other value if you - // want to use the - // ``weighted_kelly_indicator'' - // refinement criterion. - SmartPointer > kelly_weight; - - // Finally, we have a - // variable that denotes the - // maximum number of degrees - // of freedom we allow for - // the (primal) - // discretization. If it is - // exceeded, we stop the - // process of solving and - // intermittend mesh - // refinement. Its default - // value is 20,000. - unsigned int max_degrees_of_freedom; - - // Finally the default - // constructor of this class: - ProblemDescription (); - }; - - // The driver framework class - // only has one method which - // calls solver and mesh - // refinement intermittently, and - // does some other small tasks in - // between. Since it does not - // need data besides the - // parameters given to it, we - // make it static: - static void run (const ProblemDescription &descriptor); -}; - - -// As for the implementation, first -// the constructor of the parameter -// object, setting all values to -// their defaults: -template -Framework::ProblemDescription::ProblemDescription () - : - primal_fe_degree (1), - dual_fe_degree (2), - refinement_criterion (dual_weighted_error_estimator), - max_degrees_of_freedom (1000) -{} - - - -// Then the function which drives the -// whole process: -template -void Framework::run (const ProblemDescription &descriptor) -{ - // First create a triangulation - // from the given data object, - Triangulation - triangulation (Triangulation::smoothing_on_refinement); - descriptor.data->create_coarse_grid (triangulation); - - // then a set of finite elements - // and appropriate quadrature - // formula: - const FE_Q primal_fe(descriptor.primal_fe_degree); - const FE_Q dual_fe(descriptor.dual_fe_degree); - const QGauss quadrature(descriptor.dual_fe_degree+1); - const QGauss face_quadrature(descriptor.dual_fe_degree+1); - - // Next, select one of the classes - // implementing different - // refinement criteria. - LaplaceSolver::Base *solver = 0; - switch (descriptor.refinement_criterion) - { - case ProblemDescription::dual_weighted_error_estimator: - { - solver - = new LaplaceSolver::WeightedResidual (triangulation, - primal_fe, - dual_fe, - quadrature, - face_quadrature, - descriptor.data->get_right_hand_side(), - descriptor.data->get_boundary_values(), - *descriptor.dual_functional); - break; - }; - - case ProblemDescription::global_refinement: - { - solver - = new LaplaceSolver::RefinementGlobal (triangulation, - primal_fe, - quadrature, - face_quadrature, - descriptor.data->get_right_hand_side(), - descriptor.data->get_boundary_values()); - break; - }; - - case ProblemDescription::kelly_indicator: - { - solver - = new LaplaceSolver::RefinementKelly (triangulation, - primal_fe, - quadrature, - face_quadrature, - descriptor.data->get_right_hand_side(), - descriptor.data->get_boundary_values()); - break; - }; - - case ProblemDescription::weighted_kelly_indicator: - { - solver - = new LaplaceSolver::RefinementWeightedKelly (triangulation, - primal_fe, - quadrature, - face_quadrature, - descriptor.data->get_right_hand_side(), - descriptor.data->get_boundary_values(), - *descriptor.kelly_weight); - break; - }; - - default: - AssertThrow (false, ExcInternalError()); - }; - - // Now that all objects are in - // place, run the main loop. The - // stopping criterion is - // implemented at the bottom of the - // loop. - // - // In the loop, first set the new - // cycle number, then solve the - // problem, output its solution(s), - // apply the evaluation objects to - // it, then decide whether we want - // to refine the mesh further and - // solve again on this mesh, or - // jump out of the loop. - for (step=0; step<=n_steps; ++step) - { - std::cout << "Refinement cycle: " << step - << std::endl; - - solver->set_refinement_cycle (step); - solver->solve_problem (); - solver->output_solution (); - - std::cout << " Number of degrees of freedom=" - << solver->n_dofs() << std::endl; - - for (typename EvaluatorList::const_iterator - e = descriptor.evaluator_list.begin(); - e != descriptor.evaluator_list.end(); ++e) - { - (*e)->set_refinement_cycle (step); - solver->postprocess (**e); - }; - - - if (solver->n_dofs() < descriptor.max_degrees_of_freedom) - solver->refine_grid (); - else - break; - }; - - // After the loop has run, clean up - // the screen, and delete objects - // no more needed: - std::cout << std::endl; - delete solver; - solver = 0; -} - - - - -// @sect3{The main function} - -// Here finally comes the main -// function. It drives the whole -// process by specifying a set of -// parameters to be used for the -// simulation (polynomial degrees, -// evaluation and dual functionals, -// etc), and passes them packed into -// a structure to the frame work -// class above. -int main (int argc, char **argv) -{ - // if no argument is given, then do 18 - // iterations - if (argc == 1) - n_steps = 18; - else if (argc == 2) - { - n_steps = atoi(argv[1]); - if ((n_steps==0) || (n_steps>100)) - { - std::cout << "Please call this program with an argument in the range 1..100" - << std::endl; - exit (1); - } - } - else - { - std::cout << "Please call this program with a single argument in the range 1..100" - << std::endl; - exit (1); - } - - - - - deallog.depth_console (0); - try - { - // Describe the problem we want - // to solve here by passing a - // descriptor object to the - // function doing the rest of - // the work: - const unsigned int dim = 3; - Framework::ProblemDescription descriptor; - - // First set the refinement - // criterion we wish to use: - descriptor.refinement_criterion - = Framework::ProblemDescription::dual_weighted_error_estimator; - // Here, we could as well have - // used ``global_refinement'' - // or - // ``weighted_kelly_indicator''. Note - // that the information given - // about dual finite elements, - // dual functional, etc is only - // important for the given - // choice of refinement - // criterion, and is ignored - // otherwise. - - // Then set the polynomial - // degrees of primal and dual - // problem. We choose here - // bi-linear and bi-quadratic - // ones: - descriptor.primal_fe_degree = 1; - descriptor.dual_fe_degree = 2; - - // Then set the description of - // the test case, i.e. domain, - // boundary values, and right - // hand side. These are - // prepackaged in classes. We - // take here the description of - // ``Exercise_2_3'', but you - // can also use - // ``CurvedRidges'': - descriptor.data = new Data::SetUp,dim> (); - - // Next set first a dual - // functional, then a list of - // evaluation objects. We - // choose as default the - // evaluation of the - // value at an - // evaluation point, - // represented by the classes - // ``PointValueEvaluation'' - // in the namespaces of - // evaluation and dual - // functional classes. You can - // also set the - // ``PointXDerivativeEvaluation'' - // classes for the x-derivative - // instead of the value - // at the evaluation point. - // - // Note that dual functional - // and evaluation objects - // should match. However, you - // can give as many evaluation - // functionals as you want, so - // you can have both point - // value and derivative - // evaluated after each step. - // One such additional - // evaluation is to output the - // grid in each step. - const Point evaluation_point (0., 0., 0.); - descriptor.dual_functional - = new DualFunctional::PointValueEvaluation (evaluation_point); - - Evaluation::PointValueEvaluation - postprocessor1 (evaluation_point); - Evaluation::GridOutput - postprocessor2 ("grid"); - - descriptor.evaluator_list.push_back (&postprocessor1); - descriptor.evaluator_list.push_back (&postprocessor2); - - // Set the maximal number of - // degrees of freedom after - // which we want the program to - // stop refining the mesh - // further: -#if defined(SPEC_CPU) - // raise from 20000 to 30000. (jfk p6f) - descriptor.max_degrees_of_freedom = 30000; -#else - descriptor.max_degrees_of_freedom = 20000; -#endif - - // Finally pass the descriptor - // object to a function that - // runs the entire solution - // with it: - Framework::run (descriptor); - } - - // Catch exceptions to give - // information about things that - // failed: - catch (std::exception &exc) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - }; - - return 0; -} -- 2.39.5