// $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
/**
- * 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 <int dim>
+ class SlitSingularityFunction : public Function<dim>
{
public:
/**
* The value at a single point.
*/
- virtual double value (const Point<2> &p,
+ virtual double value (const Point<dim> &p,
const unsigned int component = 0) const;
/**
* Values at multiple points.
*/
- virtual void value_list (const std::vector<Point<2> > &points,
+ virtual void value_list (const std::vector<Point<dim> > &points,
std::vector<double> &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<dim> &p,
const unsigned int component = 0) const;
/**
* Gradients at multiple points.
*/
- virtual void gradient_list (const std::vector<Point<2> > &points,
- std::vector<Tensor<1,2> > &gradients,
+ virtual void gradient_list (const std::vector<Point<dim> > &points,
+ std::vector<Tensor<1,dim> > &gradients,
const unsigned int component = 0) const;
/**
* Laplacian at a single point.
*/
- virtual double laplacian (const Point<2> &p,
+ virtual double laplacian (const Point<dim> &p,
const unsigned int component = 0) const;
/**
* Laplacian at multiple points.
*/
- virtual void laplacian_list (const std::vector<Point<2> > &points,
+ virtual void laplacian_list (const std::vector<Point<dim> > &points,
std::vector<double> &values,
const unsigned int component = 0) const;
};
// $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
//////////////////////////////////////////////////////////////////////
-
+ template <int dim>
double
- SlitSingularityFunction::value (const Point<2> &p,
- const unsigned int) const
+ SlitSingularityFunction<dim>::value (
+ const Point<dim> &p,
+ const unsigned int) const
{
double x = p(0);
double y = p(1);
}
+ template <int dim>
void
- SlitSingularityFunction::value_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
- const unsigned int) const
+ SlitSingularityFunction<dim>::value_list (
+ const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
{
Assert (values.size() == points.size(),
ExcDimensionMismatch(values.size(), points.size()));
}
+ template <int dim>
double
- SlitSingularityFunction::laplacian (const Point<2> &,
+ SlitSingularityFunction<dim>::laplacian (const Point<dim> &,
const unsigned int) const
{
return 0.;
}
+ template <int dim>
void
- SlitSingularityFunction::laplacian_list (const std::vector<Point<2> > &points,
- std::vector<double> &values,
- const unsigned int) const
+ SlitSingularityFunction<dim>::laplacian_list (
+ const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int) const
{
Assert (values.size() == points.size(),
ExcDimensionMismatch(values.size(), points.size()));
}
- Tensor<1,2>
- SlitSingularityFunction::gradient (const Point<2> &p,
+ template <int dim>
+ Tensor<1,dim>
+ SlitSingularityFunction<dim>::gradient (const Point<dim> &p,
const unsigned int) const
{
double x = p(0);
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 <int dim>
void
- SlitSingularityFunction::gradient_list (const std::vector<Point<2> > &points,
- std::vector<Tensor<1,2> > &gradients,
+ SlitSingularityFunction<dim>::gradient_list (const std::vector<Point<dim> > &points,
+ std::vector<Tensor<1,dim> > &gradients,
const unsigned int) const
{
Assert (gradients.size() == points.size(),
for (unsigned int i=0;i<points.size();++i)
{
- const Point<2>& p = points[i];
+ const Point<dim>& p = points[i];
double x = p(0);
double y = p(1);
double phi = std::atan2(x,y)+M_PI;
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<dim;++d)
+ gradients[i][d] = 0.;
}
}
template class FourierSineSum<1>;
template class FourierSineSum<2>;
template class FourierSineSum<3>;
+ template class SlitSingularityFunction<2>;
+ template class SlitSingularityFunction<3>;
}
* the square at
* <tt>x=y=(left+right)/2</tt>.
*
- * The triangulation needs to be
- * void upon calling this
- * function.
- *
+ * In 3d, the 2d domain is just
+ * extended in the
+ * <i>z</i>-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
// Implementation for 3D only
template <int dim>
-void GridGenerator::hyper_cube_slit (Triangulation<dim> &,
- const double,
- const double,
- const bool)
+void GridGenerator::hyper_cube_slit (Triangulation<dim>& 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<dim> vertices[20] = {
+ Point<dim>(left, left , -len/2.),
+ Point<dim>(rl2, left , -len/2.),
+ Point<dim>(rl2, rl2 , -len/2.),
+ Point<dim>(left, rl2 , -len/2.),
+ Point<dim>(right,left , -len/2.),
+ Point<dim>(right,rl2 , -len/2.),
+ Point<dim>(rl2, right, -len/2.),
+ Point<dim>(left, right, -len/2.),
+ Point<dim>(right,right, -len/2.),
+ Point<dim>(rl2, left , -len/2.),
+ Point<dim>(left, left , len/2.),
+ Point<dim>(rl2, left , len/2.),
+ Point<dim>(rl2, rl2 , len/2.),
+ Point<dim>(left, rl2 , len/2.),
+ Point<dim>(right,left , len/2.),
+ Point<dim>(right,rl2 , len/2.),
+ Point<dim>(rl2, right, len/2.),
+ Point<dim>(left, right, len/2.),
+ Point<dim>(right,right, len/2.),
+ Point<dim>(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<CellData<dim> > cells (4, CellData<dim>());
+ 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<Point<dim> >(&vertices[0], &vertices[20]),
+ cells,
+ SubCellData()); // no boundary information
+
+ if (colorize)
+ {
+ Assert(false, ExcNotImplemented());
+ typename Triangulation<dim>::cell_iterator cell = tria.begin();
+ cell->face(1)->set_boundary_indicator(1);
+ ++cell;
+ cell->face(3)->set_boundary_indicator(2);
+ }
}