From 3180b833c96e443c688455930f2ab42f95329f0a Mon Sep 17 00:00:00 2001 From: guido Date: Wed, 7 Jun 2006 14:00:00 +0000 Subject: [PATCH] extend slit to 3D git-svn-id: https://svn.dealii.org/trunk@13190 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/base/include/base/function_lib.h | 23 ++++---- deal.II/base/source/function_lib.cc | 47 +++++++++------ deal.II/deal.II/include/grid/grid_generator.h | 13 +++-- deal.II/deal.II/source/grid/grid_generator.cc | 58 +++++++++++++++++-- 4 files changed, 101 insertions(+), 40 deletions(-) diff --git a/deal.II/base/include/base/function_lib.h b/deal.II/base/include/base/function_lib.h index c8ce204000..94f5aba243 100644 --- a/deal.II/base/include/base/function_lib.h +++ b/deal.II/base/include/base/function_lib.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -384,50 +384,51 @@ namespace Functions /** - * Singularity on the slit domain in 2D. + * Singularity on the slit domain in 2D and 3D. * * @ingroup functions - * @author Guido Kanschat, 1999 + * @author Guido Kanschat, 1999, 2006 */ - class SlitSingularityFunction : public Function<2> + template + class SlitSingularityFunction : public Function { public: /** * The value at a single point. */ - virtual double value (const Point<2> &p, + virtual double value (const Point &p, const unsigned int component = 0) const; /** * Values at multiple points. */ - virtual void value_list (const std::vector > &points, + virtual void value_list (const std::vector > &points, std::vector &values, const unsigned int component = 0) const; /** * Gradient at a single point. */ - virtual Tensor<1,2> gradient (const Point<2> &p, + virtual Tensor<1,dim> gradient (const Point &p, const unsigned int component = 0) const; /** * Gradients at multiple points. */ - virtual void gradient_list (const std::vector > &points, - std::vector > &gradients, + virtual void gradient_list (const std::vector > &points, + std::vector > &gradients, const unsigned int component = 0) const; /** * Laplacian at a single point. */ - virtual double laplacian (const Point<2> &p, + virtual double laplacian (const Point &p, const unsigned int component = 0) const; /** * Laplacian at multiple points. */ - virtual void laplacian_list (const std::vector > &points, + virtual void laplacian_list (const std::vector > &points, std::vector &values, const unsigned int component = 0) const; }; diff --git a/deal.II/base/source/function_lib.cc b/deal.II/base/source/function_lib.cc index 48fbd9fcd4..cbecc2968d 100644 --- a/deal.II/base/source/function_lib.cc +++ b/deal.II/base/source/function_lib.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -891,10 +891,11 @@ namespace Functions ////////////////////////////////////////////////////////////////////// - + template double - SlitSingularityFunction::value (const Point<2> &p, - const unsigned int) const + SlitSingularityFunction::value ( + const Point &p, + const unsigned int) const { double x = p(0); double y = p(1); @@ -906,10 +907,12 @@ namespace Functions } + template void - SlitSingularityFunction::value_list (const std::vector > &points, - std::vector &values, - const unsigned int) const + SlitSingularityFunction::value_list ( + const std::vector > &points, + std::vector &values, + const unsigned int) const { Assert (values.size() == points.size(), ExcDimensionMismatch(values.size(), points.size())); @@ -927,18 +930,21 @@ namespace Functions } + template double - SlitSingularityFunction::laplacian (const Point<2> &, + SlitSingularityFunction::laplacian (const Point &, const unsigned int) const { return 0.; } + template void - SlitSingularityFunction::laplacian_list (const std::vector > &points, - std::vector &values, - const unsigned int) const + SlitSingularityFunction::laplacian_list ( + const std::vector > &points, + std::vector &values, + const unsigned int) const { Assert (values.size() == points.size(), ExcDimensionMismatch(values.size(), points.size())); @@ -948,8 +954,9 @@ namespace Functions } - Tensor<1,2> - SlitSingularityFunction::gradient (const Point<2> &p, + template + Tensor<1,dim> + SlitSingularityFunction::gradient (const Point &p, const unsigned int) const { double x = p(0); @@ -957,17 +964,17 @@ namespace Functions double phi = std::atan2(x,y)+M_PI; double r64 = std::pow(x*x+y*y,3./4.); - - Tensor<1,2> result; + Tensor<1,dim> result; result[0] = 1./2.*(std::sin(1./2.*phi)*x + std::cos(1./2.*phi)*y)/r64; result[1] = 1./2.*(std::sin(1./2.*phi)*y - std::cos(1./2.*phi)*x)/r64; return result; } + template void - SlitSingularityFunction::gradient_list (const std::vector > &points, - std::vector > &gradients, + SlitSingularityFunction::gradient_list (const std::vector > &points, + std::vector > &gradients, const unsigned int) const { Assert (gradients.size() == points.size(), @@ -975,7 +982,7 @@ namespace Functions for (unsigned int i=0;i& p = points[i]; + const Point& p = points[i]; double x = p(0); double y = p(1); double phi = std::atan2(x,y)+M_PI; @@ -983,6 +990,8 @@ namespace Functions gradients[i][0] = 1./2.*(std::sin(1./2.*phi)*x + std::cos(1./2.*phi)*y)/r64; gradients[i][1] = 1./2.*(std::sin(1./2.*phi)*y - std::cos(1./2.*phi)*x)/r64; + for (unsigned int d=2;d; template class FourierSineSum<2>; template class FourierSineSum<3>; + template class SlitSingularityFunction<2>; + template class SlitSingularityFunction<3>; } diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index c4ea501272..f5caf5fec2 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -421,16 +421,17 @@ class GridGenerator * the square at * x=y=(left+right)/2. * - * The triangulation needs to be - * void upon calling this - * function. - * + * In 3d, the 2d domain is just + * extended in the + * z-direction, such that + * a plane cuts the lower half of + * a rectangle in two. + * This function is declared to * exist for triangulations of * all space dimensions, but * throws an error if called in - * 1d. It is also presently not - * implemented in 3d. + * 1d. * * @note The triangulation needs to be * void upon calling this diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index f2e46282be..42dbef8297 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -1037,12 +1037,60 @@ GridGenerator::half_hyper_shell (Triangulation &tria, // Implementation for 3D only template -void GridGenerator::hyper_cube_slit (Triangulation &, - const double, - const double, - const bool) +void GridGenerator::hyper_cube_slit (Triangulation& tria, + const double left, + const double right, + const bool colorize) { - Assert (false, ExcNotImplemented()); + const double rl2=(right+left)/2; + const double len = (right-left)/2.; + + const Point vertices[20] = { + Point(left, left , -len/2.), + Point(rl2, left , -len/2.), + Point(rl2, rl2 , -len/2.), + Point(left, rl2 , -len/2.), + Point(right,left , -len/2.), + Point(right,rl2 , -len/2.), + Point(rl2, right, -len/2.), + Point(left, right, -len/2.), + Point(right,right, -len/2.), + Point(rl2, left , -len/2.), + Point(left, left , len/2.), + Point(rl2, left , len/2.), + Point(rl2, rl2 , len/2.), + Point(left, rl2 , len/2.), + Point(right,left , len/2.), + Point(right,rl2 , len/2.), + Point(rl2, right, len/2.), + Point(left, right, len/2.), + Point(right,right, len/2.), + Point(rl2, left , len/2.) + }; + const int cell_vertices[4][8] = { { 0,1,3,2, 10, 11, 13, 12 }, + { 9,4,2,5, 19,14, 12, 15 }, + { 3,2,7,6,13,12,17,16 }, + { 2,5,6,8,12,15,16,18 } }; + std::vector > cells (4, CellData()); + for (unsigned int i=0; i<4; ++i) + { + for (unsigned int j=0; j<8; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + }; + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[20]), + cells, + SubCellData()); // no boundary information + + if (colorize) + { + Assert(false, ExcNotImplemented()); + typename Triangulation::cell_iterator cell = tria.begin(); + cell->face(1)->set_boundary_indicator(1); + ++cell; + cell->face(3)->set_boundary_indicator(2); + } } -- 2.39.5