<h3>Specific improvements</h3>
<ol>
+<li> New: The GridGenerator::torus function and TorusBoundary class can
+create and describe the surface of a torus.
+<br>
+(Daniel Castanon Quiroz, 2011/09/08)
+
<li> Fixed: FEFieldFunction class was not thread-safe because it keeps a
cache on the side that was invalidated when different
threads kept pouncing on it. This is now fixed.
<br>
(Patrick Sodré, 2011/09/07)
-
<li> New: There is now a function GridTools::volume computing the volume
of a triangulation.
<br>
*
*
* @ingroup grid
- * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, Stefan
- * Nauber, Joerg Weimar, Yaqi Wang, Luca Heltai, 1998, 1999, 2000, 2001, 2002,
- * 2003, 2006, 2007, 2008, 2009, 2010, 2011.
*/
class GridGenerator
{
* computed adaptively such that
* the resulting elements have
* the least aspect ratio.
- *
+ *
* If colorize is set to true, the
* inner, outer, left, and right
* boundary get indicator 0, 1, 2,
- * and 3, respectively. Otherwise
+ * and 3, respectively. Otherwise
* all indicators are set to 0.
*
* @note The triangulation needs to be
const unsigned int n_cells = 0,
const bool colorize = false);
-
+
/**
* Produce a quarter hyper-shell,
* i.e. the space between two
const double outer_radius,
const unsigned int n_cells = 0,
const bool colorize = false);
-
+
/**
* Produce a domain that is the space
* between two cylinders in 3d, with
const unsigned int n_radial_cells = 0,
const unsigned int n_axial_cells = 0);
+
+
+ /**
+ * Produce the surface meshing of the
+ * torus. The axis of the torus is the
+ * $y$-axis while the plane of the torus
+ * is the $x$-$z$ plane. The boundary of
+ * this object can be described by the
+ * TorusBoundary class.
+ *
+ * @param R The radius of the circle,
+ * which forms the middle line of the
+ * torus containing the loop of
+ * cells. Must be greater than @p r.
+ *
+ * @param r The inner radius of the
+ * torus.
+ */
+
+ static void torus (Triangulation<2,3>& tria,
+ const double R,
+ const double r);
+
+
/**
* This class produces a square
* on the <i>xy</i>-plane with a
* the third parameter. Previous
* content of @p result will be
* deleted.
- *
- * This function is most often used
+ *
+ * This function is most often used
* to compose meshes for more
* complicated geometries if the
* geometry can be composed of
* exist to generate coarse meshes.
* For example, the channel mesh used
* in step-35 could in principle be
- * created using a mesh created by the
+ * created using a mesh created by the
* GridGenerator::hyper_cube_with_cylindrical_hole
* function and several rectangles,
* and merging them using the current
* individual mesh building blocks are
* GridTools::transform, GridTools::rotate,
* and GridTools::scale).
- *
+ *
* @note The two input triangulations
* must be coarse meshes that have
- * no refined cells.
- *
+ * no refined cells.
+ *
* @note The function copies the material ids
* of the cells of the two input
* triangulations into the output
* then you will currently have to set
* boundary indicators again by hand
* in the output triangulation.
- *
+ *
* @note For a related operation
* on refined meshes when both
* meshes are derived from the
* this class and the documentation of the
* base class.
*/
- virtual
+ virtual
Point<spacedim>
get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
* this class and the documentation of the
* base class.
*/
- virtual
+ virtual
Point<spacedim>
get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
* Calls
* @p get_intermediate_points_between_points.
*/
- virtual
+ virtual
void
get_intermediate_points_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line,
std::vector<Point<spacedim> > &points) const;
* Only implemented for <tt>dim=3</tt>
* and for <tt>points.size()==1</tt>.
*/
- virtual
+ virtual
void
get_intermediate_points_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad,
std::vector<Point<spacedim> > &points) const;
* and the documentation of the
* base class.
*/
- virtual
+ virtual
void
get_normals_at_vertices (const typename Triangulation<dim,spacedim>::face_iterator &face,
typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const;
/**
* Return the radius of the ball.
*/
- double
+ double
get_radius () const;
/**
};
+/**
+ * Class describing the boundary of the torus. The axis of the torus is the
+ * $y$-axis while the plane of the torus is the $x$-$z$ plane. A torus of this
+ * kind can be generated by GridGenerator::torus.
+ *
+ * This class is only implemented for
+ * <tt>dim=2</tt>,<tt>spacedim=3</tt>, that is, just the surface.
+ */
+template <int dim, int spacedim>
+class TorusBoundary : public Boundary<dim,spacedim>
+{
+ public:
+ /**
+ * Constructor.<tt>R</tt> has to be
+ * greater than <tt>r</tt>.
+ */
+ TorusBoundary (const double R, const double r);
+
+//Boundary Refinenment Functions
+ /**
+ * Construct a new point on a line.
+ */
+ virtual Point<spacedim>
+ get_new_point_on_line (const typename Triangulation<dim,spacedim>::line_iterator &line) const;
+
+ /**
+ * Construct a new point on a quad.
+ */
+ virtual Point<spacedim>
+ get_new_point_on_quad (const typename Triangulation<dim,spacedim>::quad_iterator &quad) const;
+
+ /**
+ * Construct a new points on a line.
+ */
+ virtual void get_intermediate_points_on_line (
+ const typename Triangulation< dim, spacedim >::line_iterator & line,
+ std::vector< Point< spacedim > > & points) const ;
+
+ /**
+ * Construct a new points on a quad.
+ */
+ virtual void get_intermediate_points_on_quad (
+ const typename Triangulation< dim, spacedim >::quad_iterator & quad,
+ std::vector< Point< spacedim > > & points ) const ;
+
+ /**
+ * Get the normal from cartesian
+ * coordinates. This normal does not have
+ * unit length.
+ */
+ virtual void get_normals_at_vertices (
+ const typename Triangulation< dim, spacedim >::face_iterator &face,
+ typename Boundary<dim,spacedim>::FaceVertexNormals &face_vertex_normals) const ;
+
+ private:
+ //Handy functions
+ /**
+ * Function that corrects the value and
+ * sign of angle, that is, given
+ * <tt>angle=tan(abs(y/x))<tt>; if <tt>
+ * (y > 0) && (x < 0) <tt> then
+ * <tt>correct_angle = Pi - angle<tt>,
+ * etc.
+ */
+
+ double get_correct_angle(const double angle,const double x,const double y) const;
+
+ /**
+ * Get the cartesian coordinates of the
+ * Torus, i.e., from <tt>(theta,phi)<tt>
+ * to <tt>(x,y,z)<tt>.
+ */
+ Point<spacedim> get_real_coord(const Point<dim>& surfP) const ;
+
+ /**
+ * Get the surface coordinates of the
+ * Torus, i.e., from <tt>(x,y,z)<tt> to
+ * <tt>(theta,phi)<tt>.
+ */
+ Point<dim> get_surf_coord(const Point<spacedim>& p) const ;
+
+ /**
+ * Get the normal from surface
+ * coordinates. This normal does not have
+ * unit length.
+ */
+ Point<spacedim> get_surf_norm_from_sp(const Point<dim>& surfP) const ;
+
+ /**
+ * Get the normal from cartesian
+ * coordinates. This normal does not have
+ * unit length.
+ */
+ Point<spacedim> get_surf_norm(const Point<spacedim>& p) const ;
+
+ /**
+ * Inner and outer radii of the shell.
+ */
+ const double R;
+ const double r;
+};
+
+
/* -------------- declaration of explicit specializations ------------- */
get_normals_at_vertices (const Triangulation<1,1>::face_iterator &,
Boundary<1,1>::FaceVertexNormals &) const;
+
#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE
+void
+GridGenerator::torus (Triangulation<2,3>& tria,
+ const double R,
+ const double r)
+{
+ Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
+
+ const unsigned int dim=2;
+ const unsigned int spacedim=3;
+ std::vector<Point<spacedim> > vertices (16);
+
+ vertices[0]=Point<spacedim>(R-r,0,0);
+ vertices[1]=Point<spacedim>(R,-r,0);
+ vertices[2]=Point<spacedim>(R+r,0,0);
+ vertices[3]=Point<spacedim>(R, r,0);
+ vertices[4]=Point<spacedim>(0,0,R-r);
+ vertices[5]=Point<spacedim>(0,-r,R);
+ vertices[6]=Point<spacedim>(0,0,R+r);
+ vertices[7]=Point<spacedim>(0,r,R);
+ vertices[8]=Point<spacedim>(-(R-r),0,0);
+ vertices[9]=Point<spacedim>(-R,-r,0);
+ vertices[10]=Point<spacedim>(-(R+r),0,0);
+ vertices[11]=Point<spacedim>(-R, r,0);
+ vertices[12]=Point<spacedim>(0,0,-(R-r));
+ vertices[13]=Point<spacedim>(0,-r,-R);
+ vertices[14]=Point<spacedim>(0,0,-(R+r));
+ vertices[15]=Point<spacedim>(0,r,-R);
+
+ std::vector<CellData<dim> > cells (16);
+ //Right Hand Orientation
+ cells[0].vertices[0] = 0;
+ cells[0].vertices[1] = 4;
+ cells[0].vertices[2] = 7;
+ cells[0].vertices[3] = 3;
+ cells[0].material_id = 0;
+
+ cells[1].vertices[0] = 1;
+ cells[1].vertices[1] = 5;
+ cells[1].vertices[2] = 4;
+ cells[1].vertices[3] = 0;
+ cells[1].material_id = 0;
+
+ cells[2].vertices[0] = 2;
+ cells[2].vertices[1] = 6;
+ cells[2].vertices[2] = 5;
+ cells[2].vertices[3] = 1;
+ cells[2].material_id = 0;
+
+ cells[3].vertices[0] = 3;
+ cells[3].vertices[1] = 7;
+ cells[3].vertices[2] = 6;
+ cells[3].vertices[3] = 2;
+ cells[3].material_id = 0;
+
+ cells[4].vertices[0] = 4;
+ cells[4].vertices[1] = 8;
+ cells[4].vertices[2] = 11;
+ cells[4].vertices[3] = 7;
+ cells[4].material_id = 0;
+
+ cells[5].vertices[0] = 5;
+ cells[5].vertices[1] = 9;
+ cells[5].vertices[2] = 8;
+ cells[5].vertices[3] = 4;
+ cells[5].material_id = 0;
+
+ cells[6].vertices[0] = 6;
+ cells[6].vertices[1] = 10;
+ cells[6].vertices[2] = 9;
+ cells[6].vertices[3] = 5;
+ cells[6].material_id = 0;
+
+ cells[7].vertices[0] = 7;
+ cells[7].vertices[1] = 11;
+ cells[7].vertices[2] = 10;
+ cells[7].vertices[3] = 6;
+ cells[7].material_id = 0;
+
+ cells[8].vertices[0] = 8;
+ cells[8].vertices[1] = 12;
+ cells[8].vertices[2] = 15;
+ cells[8].vertices[3] = 11;
+ cells[8].material_id = 0;
+
+ cells[9].vertices[0] = 9;
+ cells[9].vertices[1] = 13;
+ cells[9].vertices[2] = 12;
+ cells[9].vertices[3] = 8;
+ cells[9].material_id = 0;
+
+ cells[10].vertices[0] = 10;
+ cells[10].vertices[1] = 14;
+ cells[10].vertices[2] = 13;
+ cells[10].vertices[3] = 9;
+ cells[10].material_id = 0;
+
+ cells[11].vertices[0] = 11;
+ cells[11].vertices[1] = 15;
+ cells[11].vertices[2] = 14;
+ cells[11].vertices[3] = 10;
+ cells[11].material_id = 0;
+
+ cells[12].vertices[0] = 12;
+ cells[12].vertices[1] = 0;
+ cells[12].vertices[2] = 3;
+ cells[12].vertices[3] = 15;
+ cells[12].material_id = 0;
+
+ cells[13].vertices[0] = 13;
+ cells[13].vertices[1] = 1;
+ cells[13].vertices[2] = 0;
+ cells[13].vertices[3] = 12;
+ cells[13].material_id = 0;
+
+ cells[14].vertices[0] = 14;
+ cells[14].vertices[1] = 2;
+ cells[14].vertices[2] = 1;
+ cells[14].vertices[3] = 13;
+ cells[14].material_id = 0;
+
+ cells[15].vertices[0] = 15;
+ cells[15].vertices[1] = 3;
+ cells[15].vertices[2] = 2;
+ cells[15].vertices[3] = 14;
+ cells[15].material_id = 0;
+
+ // Must call this to be able to create a
+ // correct triangulation in dealii, read
+ // GridReordering<> doc
+ GridReordering<dim,spacedim>::reorder_cells (cells);
+ tria.create_triangulation_compatibility (vertices, cells, SubCellData());
+}
+
+
+
// Implementation for 2D only
template<>
};
tria.create_triangulation (vertices, cells, SubCellData());
-
+
if (colorize)
{
Triangulation<2>::cell_iterator cell = tria.begin();
cell->face(2)->set_boundary_indicator(1);
}
tria.begin()->face(0)->set_boundary_indicator(3);
-
+
tria.last()->face(1)->set_boundary_indicator(2);
}
}
cells[i].material_id = 0;
};
- tria.create_triangulation (vertices, cells, SubCellData());
-
+ tria.create_triangulation (vertices, cells, SubCellData());
+
if (colorize)
{
Triangulation<2>::cell_iterator cell = tria.begin();
cell->face(2)->set_boundary_indicator(1);
}
tria.begin()->face(0)->set_boundary_indicator(3);
-
+
tria.last()->face(1)->set_boundary_indicator(2);
}
}
ExcMessage ("The input triangulations must be coarse meshes."));
Assert (triangulation_2.n_levels() == 1,
ExcMessage ("The input triangulations must be coarse meshes."));
-
+
// get the union of the set of vertices
std::vector<Point<spacedim> > vertices = triangulation_1.get_vertices();
vertices.insert (vertices.end(),
triangulation_2.get_vertices().begin(),
triangulation_2.get_vertices().end());
-
+
// now form the union of the set of cells
std::vector<CellData<dim> > cells;
cells.reserve (triangulation_1.n_cells() + triangulation_2.n_cells());
this_cell.material_id = cell->material_id();
cells.push_back (this_cell);
}
-
- // now do the same for the other other mesh. note that we have to
+
+ // now do the same for the other other mesh. note that we have to
// translate the vertex indices
for (typename Triangulation<dim,spacedim>::cell_iterator
cell = triangulation_2.begin(); cell != triangulation_2.end(); ++cell)
this_cell.material_id = cell->material_id();
cells.push_back (this_cell);
}
-
+
// throw out duplicated vertices from the two meshes
// and create the triangulation
SubCellData subcell_data;
+
+template <int dim, int spacedim>
+TorusBoundary<dim,spacedim>::TorusBoundary (const double R__,
+ const double r__)
+ :
+ R(R__),
+ r(r__)
+{
+ Assert (false, ExcNotImplemented());
+}
+
+
+
+template <>
+TorusBoundary<2,3>::TorusBoundary (const double R__,
+ const double r__)
+ :
+ R(R__),
+ r(r__)
+{
+ Assert (R>r, ExcMessage("Outer radius must be greater than inner radius."));
+}
+
+
+
+template <int dim, int spacedim>
+double
+TorusBoundary<dim,spacedim>::get_correct_angle(const double angle,
+ const double x,
+ const double y) const
+{
+ if(y>=0)
+ {
+ if (x >=0)
+ return angle;
+
+ return numbers::PI-angle;
+ }
+
+ if (x <=0)
+ return numbers::PI+angle;
+
+ return 2.0*numbers::PI-angle;
+}
+
+
+
+template <>
+Point<3>
+TorusBoundary<2,3>::get_real_coord (const Point<2>& surfP) const
+{
+ const double theta=surfP(0);
+ const double phi=surfP(1);
+
+ return Point<3> ((R+r*std::cos(phi))*std::cos(theta),
+ r*std::sin(phi),
+ (R+r*std::cos(phi))*std::sin(theta));
+}
+
+
+
+template <>
+Point<2>
+TorusBoundary<2,3>::get_surf_coord(const Point<3>& p) const
+{
+ const double phi=std::asin(std::abs(p(1))/r);
+ const double Rr_2=p(0)*p(0)+p(2)*p(2);
+
+ Point<2> surfP;
+ surfP(1)=get_correct_angle(phi,Rr_2-R*R,p(1));//phi
+
+ if(std::abs(p(0))<1.E-5)
+ {
+ if(p(2)>=0)
+ surfP(0) = numbers::PI*0.5;
+ else
+ surfP(0) = -numbers::PI*0.5;
+ }
+ else
+ {
+ const double theta = std::atan(std::abs(p(2)/p(0)));
+ surfP(0)= get_correct_angle(theta,p(0),p(2));
+ }
+
+ return surfP;
+}
+
+
+
+template <>
+Point<3>
+TorusBoundary<2,3>::get_new_point_on_line (const typename Triangulation<2,3>::line_iterator &line) const
+{
+ //Just get the average
+ Point<2> p0=get_surf_coord(line->vertex(0));
+ Point<2> p1=get_surf_coord(line->vertex(1));
+
+ Point<2> middle(0,0);
+
+ //Take care for periodic conditions,
+ //For instance phi0= 0, phi1= 3/2*Pi middle has to be 7/4*Pi not 3/4*Pi
+ //This also works for -Pi/2 + Pi, middle is 5/4*Pi
+ for(unsigned i=0; i<2;i++)
+ if(std::abs(p0(i)-p1(i))> numbers::PI)
+ middle(i)=2*numbers::PI;
+
+ middle+= p0 + p1;
+ middle*=0.5;
+
+ Point<3> midReal=get_real_coord(middle);
+ return midReal;
+}
+
+
+
+template <>
+Point<3>
+TorusBoundary<2,3>::get_new_point_on_quad (const typename Triangulation<2,3>::quad_iterator &quad) const
+{
+ //Just get the average
+ Point<2> p[4];
+
+ for(unsigned i=0;i<4;i++)
+ p[i]=get_surf_coord(quad->vertex(i));
+
+ Point<2> middle(0,0);
+
+ //Take care for periodic conditions, see get_new_point_on_line() above
+ //For instance phi0= 0, phi1= 3/2*Pi middle has to be 7/4*Pi not 3/4*Pi
+ //This also works for -Pi/2 + Pi + Pi- Pi/2, middle is 5/4*Pi
+ for(unsigned i=0;i<2;i++)
+ for(unsigned j=1;j<4;j++){
+ if(std::abs(p[0](i)-p[j](i))> numbers::PI){
+ middle(i)+=2*numbers::PI;
+ }
+ }
+
+ for(unsigned i=0;i<4;i++)
+ middle+=p[i];
+
+ middle*= 0.25;
+
+ return get_real_coord(middle);
+}
+
+
+
+//Normal field without unit length
+template <>
+Point<3>
+TorusBoundary<2,3>:: get_surf_norm_from_sp(const Point<2>& surfP) const
+{
+
+ Point<3> n;
+ double theta=surfP[0];
+ double phi=surfP[1];
+
+ double f=R+r*std::cos(phi);
+
+ n[0]=r*std::cos(phi)*std::cos(theta)*f;
+ n[1]=r*std::sin(phi)*f;
+ n[2]=r*std::sin(theta)*std::cos(phi)*f;
+
+ return n;
+}
+
+
+
+//Normal field without unit length
+template <>
+Point<3>
+TorusBoundary<2,3>::get_surf_norm(const Point<3>& p) const{
+
+ Point<2> surfP=get_surf_coord(p);
+ return get_surf_norm_from_sp(surfP);
+
+}
+
+
+
+template<>
+void
+TorusBoundary<2,3>::
+get_intermediate_points_on_line (const typename Triangulation<2, 3>::line_iterator & line,
+ std::vector< Point< 3 > > & points) const
+{
+ //Almost the same implementation as StraightBoundary<2,3>
+ unsigned npoints=points.size();
+ if(npoints==0) return;
+
+ Point<2> p[2];
+
+ for(unsigned i=0;i<2;i++)
+ p[i]=get_surf_coord(line->vertex(i));
+
+ unsigned offset[2];
+ offset[0]=0;
+ offset[1]=0;
+
+ //Take care for periodic conditions & negative angles,
+ //see get_new_point_on_line() above
+ //Because we dont have a symmetric interpolation (just the middle) we need to
+ //add 2*Pi to each almost zero and negative angles.
+ for(unsigned i=0;i<2;i++)
+ for(unsigned j=1;j<2;j++){
+ if(std::abs(p[0](i)-p[j](i))> numbers::PI){
+ offset[i]++;
+ break;
+ }
+ }
+
+ for(unsigned i=0;i<2;i++)
+ for(unsigned j=0;j<2;j++)
+ if(p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
+ p[j](i)+=2*numbers::PI*offset[i];
+
+
+ double dx=1.0/(npoints+1);
+ double x=dx;
+ Point<2> target;
+ for(unsigned i=0; i<npoints;i++,x+=dx){
+ target= (1-x)*p[0] + x*p[1];
+ points[i]=get_real_coord(target);
+ }
+}
+
+
+
+template<>
+void
+TorusBoundary<2,3>::
+get_intermediate_points_on_quad (const typename Triangulation< 2, 3 >::quad_iterator &quad,
+ std::vector< Point< 3 > > & points )const
+{
+ //Almost the same implementation as StraightBoundary<2,3>
+ const unsigned int n=points.size(),
+ m=static_cast<unsigned int>(std::sqrt(static_cast<double>(n)));
+ // is n a square number
+ Assert(m*m==n, ExcInternalError());
+
+ const double ds=1./(m+1);
+ double y=ds;
+
+ Point<2> p[4];
+
+ for(unsigned i=0;i<4;i++)
+ p[i]=get_surf_coord(quad->vertex(i));
+
+ Point<2> target;
+ unsigned offset[2];
+ offset[0]=0;
+ offset[1]=0;
+
+ //Take care for periodic conditions & negative angles,
+ //see get_new_point_on_line() above
+ //Because we dont have a symmetric interpolation (just the middle) we need to
+ //add 2*Pi to each almost zero and negative angles.
+ for(unsigned i=0;i<2;i++)
+ for(unsigned j=1;j<4;j++){
+ if(std::abs(p[0](i)-p[j](i))> numbers::PI){
+ offset[i]++;
+ break;
+ }
+ }
+
+ for(unsigned i=0;i<2;i++)
+ for(unsigned j=0;j<4;j++)
+ if(p[j](i)<1.E-12 ) //Take care for periodic conditions & negative angles
+ p[j](i)+=2*numbers::PI*offset[i];
+
+ for (unsigned i=0; i<m; ++i, y+=ds)
+ {
+ double x=ds;
+ for (unsigned j=0; j<m; ++j, x+=ds){
+ target=((1-x) * p[0] +
+ x * p[1]) * (1-y) +
+ ((1-x) * p[2] +
+ x * p[3]) * y;
+
+ points[i*m+j]=get_real_coord(target);
+ }
+ }
+}
+
+
+
+template<>
+void
+TorusBoundary<2,3>::
+get_normals_at_vertices (const typename Triangulation<2,3 >::face_iterator &face,
+ FaceVertexNormals &face_vertex_normals) const
+{
+ for(unsigned i=0; i<GeometryInfo<2>::vertices_per_face; i++)
+ face_vertex_normals[i]=get_surf_norm(face->vertex(i));
+}
+
+
+
// explicit instantiations
#include "tria_boundary_lib.inst"
projects / dealii-svn.git / commitdiff