include/deal.II/base/geometry_info.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: geometry_info.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__geometry_info_h
00013 #define __deal2__geometry_info_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/point.h>
00019 
00020 DEAL_II_NAMESPACE_OPEN
00021 
00022 
00023 
00037 template <int dim>
00038 struct RefinementPossibilities
00039 {
00074     enum Possibilities
00075     {
00076           no_refinement= 0,
00077 
00078           isotropic_refinement = static_cast<unsigned char>(-1)
00079     };
00080 };
00081 
00082 
00083 
00094 template <>
00095 struct RefinementPossibilities<1>
00096 {
00131     enum Possibilities
00132     {
00133           no_refinement= 0,
00134           cut_x        = 1,
00135 
00136           isotropic_refinement = cut_x
00137     };
00138 };
00139 
00140 
00141 
00153 template <>
00154 struct RefinementPossibilities<2>
00155 {
00190     enum Possibilities
00191     {
00192           no_refinement= 0,
00193           cut_x        = 1,
00194           cut_y        = 2,
00195           cut_xy       = cut_x | cut_y,
00196 
00197           isotropic_refinement = cut_xy
00198     };
00199 };
00200 
00201 
00202 
00215 template <>
00216 struct RefinementPossibilities<3>
00217 {
00252     enum Possibilities
00253     {
00254           no_refinement= 0,
00255           cut_x        = 1,
00256           cut_y        = 2,
00257           cut_xy       = cut_x | cut_y,
00258           cut_z        = 4,
00259           cut_xz       = cut_x | cut_z,
00260           cut_yz       = cut_y | cut_z,
00261           cut_xyz      = cut_x | cut_y | cut_z,
00262 
00263           isotropic_refinement = cut_xyz
00264     };
00265 };
00266 
00267 
00268 
00279 template <int dim>
00280 class RefinementCase : public RefinementPossibilities<dim>
00281 {
00282   public:
00287     RefinementCase ();
00288 
00296    RefinementCase (const typename RefinementPossibilities<dim>::Possibilities refinement_case);
00297 
00307    explicit RefinementCase (const unsigned char refinement_case);
00308 
00334     operator unsigned char () const;
00335 
00342     RefinementCase operator | (const RefinementCase &r) const;
00343 
00350     RefinementCase operator & (const RefinementCase &r) const;
00351 
00367     RefinementCase operator ~ () const;
00368 
00369 
00379     static
00380     RefinementCase cut_axis (const unsigned int i);
00381 
00387     static std::size_t memory_consumption ();
00388 
00393     template <class Archive>
00394     void serialize(Archive & ar,
00395                    const unsigned int version);
00396 
00400     DeclException1 (ExcInvalidRefinementCase,
00401                     int,
00402                     << "The refinement flags given (" << arg1 << ") contain set bits that do not "
00403                     << "make sense for the space dimension of the object to which they are applied.");
00404 
00405   private:
00412     unsigned char value : (dim > 0 ? dim : 1);
00413 };
00414 
00415 
00416 namespace internal
00417 {
00418 
00419 
00441   template <int dim>
00442   struct SubfacePossibilities
00443   {
00449       enum Possibilities
00450       {
00451             case_none = 0,
00452 
00453             case_isotropic = static_cast<unsigned char>(-1)
00454       };
00455   };
00456 
00457 
00468   template <>
00469   struct SubfacePossibilities<0>
00470   {
00478       enum Possibilities
00479       {
00480             case_none = 0,
00481 
00482             case_isotropic = case_none
00483       };
00484   };
00485 
00486 
00487 
00499   template <>
00500   struct SubfacePossibilities<1>
00501   {
00511       enum Possibilities
00512       {
00513             case_none = 0,
00514 
00515             case_isotropic = case_none
00516       };
00517   };
00518 
00519 
00520 
00533   template <>
00534   struct SubfacePossibilities<2>
00535   {
00548       enum Possibilities
00549       {
00550             case_none = 0,
00551             case_x    = 1,
00552 
00553             case_isotropic = case_x
00554       };
00555   };
00556 
00557 
00558 
00654   template <>
00655   struct SubfacePossibilities<3>
00656   {
00667       enum Possibilities
00668       {
00669             case_none  = 0,
00670             case_x     = 1,
00671             case_x1y   = 2,
00672             case_x2y   = 3,
00673             case_x1y2y = 4,
00674             case_y     = 5,
00675             case_y1x   = 6,
00676             case_y2x   = 7,
00677             case_y1x2x = 8,
00678             case_xy    = 9,
00679 
00680             case_isotropic = case_xy
00681       };
00682   };
00683 
00684 
00685 
00686 
00695   template <int dim>
00696   class SubfaceCase : public SubfacePossibilities<dim>
00697   {
00698     public:
00707       SubfaceCase (const typename SubfacePossibilities<dim>::Possibilities subface_possibility);
00708 
00734       operator unsigned char () const;
00735 
00741       static std::size_t memory_consumption ();
00742 
00746       DeclException1 (ExcInvalidSubfaceCase,
00747                       int,
00748                       << "The subface case given (" << arg1 << ") does not make sense "
00749                       << "for the space dimension of the object to which they are applied.");
00750 
00751     private:
00758       unsigned char value : (dim == 3 ? 4 : 1);
00759   };
00760 
00761 } // namespace internal
00762 
00763 
00764 
00765 template <int dim> struct GeometryInfo;
00766 
00767 
00768 
00769 
00790 template <>
00791 struct GeometryInfo<0>
00792 {
00793 
00805     static const unsigned int max_children_per_cell = 1;
00806 
00810     static const unsigned int faces_per_cell    = 0;
00811 
00823     static const unsigned int max_children_per_face = 0;
00824 
00833     static unsigned int n_children(const RefinementCase<0> &refinement_case);
00834 
00838     static const unsigned int vertices_per_cell = 1;
00839 
00849     static const unsigned int vertices_per_face = 0;
00850 
00854     static const unsigned int lines_per_face    = 0;
00855 
00859     static const unsigned int quads_per_face    = 0;
00860 
00864     static const unsigned int lines_per_cell    = 0;
00865 
00870     static const unsigned int quads_per_cell    = 0;
00871 
00876     static const unsigned int hexes_per_cell    = 0;
00877 };
00878 
00879 
00880 
00881 
00882 
01414 template <int dim>
01415 struct GeometryInfo
01416 {
01417 
01429     static const unsigned int max_children_per_cell = 1 << dim;
01430 
01434     static const unsigned int faces_per_cell = 2 * dim;
01435 
01447     static const unsigned int max_children_per_face = GeometryInfo<dim-1>::max_children_per_cell;
01448 
01452     static const unsigned int vertices_per_cell = 1 << dim;
01453 
01458     static const unsigned int vertices_per_face = GeometryInfo<dim-1>::vertices_per_cell;
01459 
01463     static const unsigned int lines_per_face
01464     = GeometryInfo<dim-1>::lines_per_cell;
01465 
01469     static const unsigned int quads_per_face
01470     = GeometryInfo<dim-1>::quads_per_cell;
01471 
01484     static const unsigned int lines_per_cell
01485     = (2*GeometryInfo<dim-1>::lines_per_cell +
01486        GeometryInfo<dim-1>::vertices_per_cell);
01487 
01498     static const unsigned int quads_per_cell
01499     = (2*GeometryInfo<dim-1>::quads_per_cell +
01500        GeometryInfo<dim-1>::lines_per_cell);
01501 
01506     static const unsigned int hexes_per_cell
01507     = (2*GeometryInfo<dim-1>::hexes_per_cell +
01508        GeometryInfo<dim-1>::quads_per_cell);
01509 
01536     static const unsigned int ucd_to_deal[vertices_per_cell];
01537 
01557     static const unsigned int dx_to_deal[vertices_per_cell];
01558 
01569     static const unsigned int vertex_to_face[vertices_per_cell][dim];
01570 
01576     static
01577     unsigned int
01578     n_children(const RefinementCase<dim> &refinement_case);
01579 
01586     static
01587     unsigned int
01588     n_subfaces(const internal::SubfaceCase<dim> &subface_case);
01589 
01605     static
01606     double
01607     subface_ratio(const internal::SubfaceCase<dim> &subface_case,
01608                   const unsigned int subface_no);
01609 
01618     static
01619     RefinementCase<dim-1>
01620     face_refinement_case (const RefinementCase<dim> &cell_refinement_case,
01621                           const unsigned int face_no,
01622                           const bool face_orientation = true,
01623                           const bool face_flip        = false,
01624                           const bool face_rotation    = false);
01625 
01634     static
01635     RefinementCase<dim>
01636     min_cell_refinement_case_for_face_refinement
01637     (const RefinementCase<dim-1> &face_refinement_case,
01638      const unsigned int face_no,
01639      const bool face_orientation = true,
01640      const bool face_flip        = false,
01641      const bool face_rotation    = false);
01642 
01650     static
01651     RefinementCase<1>
01652     line_refinement_case(const RefinementCase<dim> &cell_refinement_case,
01653                          const unsigned int line_no);
01654 
01661     static
01662     RefinementCase<dim>
01663     min_cell_refinement_case_for_line_refinement(const unsigned int line_no);
01664 
01730     static
01731     unsigned int
01732     child_cell_on_face (const RefinementCase<dim> &ref_case,
01733                         const unsigned int face,
01734                         const unsigned int subface,
01735                         const bool face_orientation = true,
01736                         const bool face_flip        = false,
01737                         const bool face_rotation    = false,
01738                         const RefinementCase<dim-1> &face_refinement_case
01739                         = RefinementCase<dim-1>::isotropic_refinement);
01740 
01762     static
01763     unsigned int
01764     line_to_cell_vertices (const unsigned int line,
01765                            const unsigned int vertex);
01766 
01799     static
01800     unsigned int
01801     face_to_cell_vertices (const unsigned int face,
01802                            const unsigned int vertex,
01803                            const bool face_orientation = true,
01804                            const bool face_flip        = false,
01805                            const bool face_rotation    = false);
01806 
01826     static
01827     unsigned int
01828     face_to_cell_lines (const unsigned int face,
01829                         const unsigned int line,
01830                         const bool face_orientation = true,
01831                         const bool face_flip        = false,
01832                         const bool face_rotation    = false);
01833 
01847     static
01848     unsigned int
01849     standard_to_real_face_vertex (const unsigned int vertex,
01850                                   const bool face_orientation = true,
01851                                   const bool face_flip        = false,
01852                                   const bool face_rotation    = false);
01853 
01867     static
01868     unsigned int
01869     real_to_standard_face_vertex (const unsigned int vertex,
01870                                   const bool face_orientation = true,
01871                                   const bool face_flip        = false,
01872                                   const bool face_rotation    = false);
01873 
01887     static
01888     unsigned int
01889     standard_to_real_face_line (const unsigned int line,
01890                                 const bool face_orientation = true,
01891                                 const bool face_flip        = false,
01892                                 const bool face_rotation    = false);
01893 
01907     static
01908     unsigned int
01909     real_to_standard_face_line (const unsigned int line,
01910                                 const bool face_orientation = true,
01911                                 const bool face_flip        = false,
01912                                 const bool face_rotation    = false);
01913 
01921     static
01922     Point<dim>
01923     unit_cell_vertex (const unsigned int vertex);
01924 
01940     static
01941     unsigned int
01942     child_cell_from_point (const Point<dim> &p);
01943 
01958     static
01959     Point<dim>
01960     cell_to_child_coordinates (const Point<dim>          &p,
01961                                const unsigned int         child_index,
01962                                const RefinementCase<dim>  refine_case
01963                                = RefinementCase<dim>::isotropic_refinement);
01964 
01973     static
01974     Point<dim>
01975     child_to_cell_coordinates (const Point<dim>          &p,
01976                                const unsigned int         child_index,
01977                                const RefinementCase<dim>  refine_case
01978                                = RefinementCase<dim>::isotropic_refinement);
01979 
01985     static
01986     bool
01987     is_inside_unit_cell (const Point<dim> &p);
01988 
02011     static
02012     bool
02013     is_inside_unit_cell (const Point<dim> &p,
02014                          const double eps);
02015 
02022     static
02023     Point<dim>
02024     project_to_unit_cell (const Point<dim> &p);
02025 
02034     static
02035     double
02036     distance_to_unit_cell (const Point<dim> &p);
02037 
02043     static
02044     double
02045     d_linear_shape_function (const Point<dim> &xi,
02046                              const unsigned int i);
02047 
02053     static
02054     Tensor<1,dim>
02055     d_linear_shape_function_gradient (const Point<dim> &xi,
02056                                       const unsigned int i);
02057 
02154     template <int spacedim>
02155     static
02156     void
02157     alternating_form_at_vertices
02158 #ifndef DEAL_II_ARRAY_ARG_BUG
02159     (const Point<spacedim> (&vertices)[vertices_per_cell],
02160      Tensor<spacedim-dim,spacedim> (&forms)[vertices_per_cell])
02161 #else
02162     (const Point<spacedim> *vertices,
02163      Tensor<spacedim-dim,spacedim> *forms)
02164 #endif
02165       ;
02166 
02184     static const unsigned int unit_normal_direction[faces_per_cell];
02185 
02217     static const int unit_normal_orientation[faces_per_cell];
02218 
02226     static const unsigned int opposite_face[faces_per_cell];
02227 
02228 
02232     DeclException1 (ExcInvalidCoordinate,
02233                     double,
02234                     << "The coordinates must satisfy 0 <= x_i <= 1, "
02235                     << "but here we have x_i=" << arg1);
02236 
02240     DeclException3 (ExcInvalidSubface,
02241                     int, int, int,
02242                     << "RefinementCase<dim> " << arg1 << ": face " << arg2
02243                     << " has no subface " << arg3);
02244 };
02245 
02246 
02247 
02248 
02249 #ifndef DOXYGEN
02250 
02251 
02252 /* -------------- declaration of explicit specializations ------------- */
02253 
02254 #ifndef DEAL_II_MEMBER_ARRAY_SPECIALIZATION_BUG
02255 template <>
02256 const unsigned int GeometryInfo<1>::unit_normal_direction[faces_per_cell];
02257 template <>
02258 const unsigned int GeometryInfo<2>::unit_normal_direction[faces_per_cell];
02259 template <>
02260 const unsigned int GeometryInfo<3>::unit_normal_direction[faces_per_cell];
02261 template <>
02262 const unsigned int GeometryInfo<4>::unit_normal_direction[faces_per_cell];
02263 
02264 template <>
02265 const int GeometryInfo<1>::unit_normal_orientation[faces_per_cell];
02266 template <>
02267 const int GeometryInfo<2>::unit_normal_orientation[faces_per_cell];
02268 template <>
02269 const int GeometryInfo<3>::unit_normal_orientation[faces_per_cell];
02270 template <>
02271 const int GeometryInfo<4>::unit_normal_orientation[faces_per_cell];
02272 
02273 template <>
02274 const unsigned int GeometryInfo<1>::opposite_face[faces_per_cell];
02275 template <>
02276 const unsigned int GeometryInfo<2>::opposite_face[faces_per_cell];
02277 template <>
02278 const unsigned int GeometryInfo<3>::opposite_face[faces_per_cell];
02279 template <>
02280 const unsigned int GeometryInfo<4>::opposite_face[faces_per_cell];
02281 #endif
02282 
02283 
02284 template <>
02285 Tensor<1,1>
02286 GeometryInfo<1>::
02287 d_linear_shape_function_gradient (const Point<1> &xi,
02288                                   const unsigned int i);
02289 template <>
02290 Tensor<1,2>
02291 GeometryInfo<2>::
02292 d_linear_shape_function_gradient (const Point<2> &xi,
02293                                   const unsigned int i);
02294 template <>
02295 Tensor<1,3>
02296 GeometryInfo<3>::
02297 d_linear_shape_function_gradient (const Point<3> &xi,
02298                                   const unsigned int i);
02299 
02300 
02301 
02302 
02303 /* -------------- inline functions ------------- */
02304 
02305 namespace internal
02306 {
02307 
02308   template <int dim>
02309   inline
02310   SubfaceCase<dim>::SubfaceCase (const typename SubfacePossibilities<dim>::Possibilities subface_possibility)
02311                   :
02312                   value (subface_possibility)
02313   {}
02314 
02315 
02316   template <int dim>
02317   inline
02318   SubfaceCase<dim>::operator unsigned char () const
02319   {
02320     return value;
02321   }
02322 
02323 
02324 } // namespace internal
02325 
02326 
02327 template <int dim>
02328 inline
02329 RefinementCase<dim>
02330 RefinementCase<dim>::cut_axis (const unsigned int)
02331 {
02332   Assert (false, ExcInternalError());
02333   return static_cast<unsigned char>(-1);
02334 }
02335 
02336 
02337 template <>
02338 inline
02339 RefinementCase<1>
02340 RefinementCase<1>::cut_axis (const unsigned int i)
02341 {
02342   const unsigned int dim = 1;
02343   Assert (i < dim, ExcIndexRange(i, 0, dim));
02344 
02345   static const RefinementCase options[dim] = { cut_x };
02346   return options[i];
02347 }
02348 
02349 
02350 
02351 template <>
02352 inline
02353 RefinementCase<2>
02354 RefinementCase<2>::cut_axis (const unsigned int i)
02355 {
02356   const unsigned int dim = 2;
02357   Assert (i < dim, ExcIndexRange(i, 0, dim));
02358 
02359   static const RefinementCase options[dim] = { cut_x, cut_y };
02360   return options[i];
02361 }
02362 
02363 
02364 
02365 template <>
02366 inline
02367 RefinementCase<3>
02368 RefinementCase<3>::cut_axis (const unsigned int i)
02369 {
02370   const unsigned int dim = 3;
02371   Assert (i < dim, ExcIndexRange(i, 0, dim));
02372 
02373   static const RefinementCase options[dim] = { cut_x, cut_y, cut_z };
02374   return options[i];
02375 }
02376 
02377 
02378 
02379 template <int dim>
02380 inline
02381 RefinementCase<dim>::RefinementCase ()
02382                 :
02383                 value (RefinementPossibilities<dim>::no_refinement)
02384 {}
02385 
02386 
02387 
02388 template <int dim>
02389 inline
02390 RefinementCase<dim>::
02391 RefinementCase (const typename RefinementPossibilities<dim>::Possibilities refinement_case)
02392                 :
02393                 value (refinement_case)
02394 {
02395                                    // check that only those bits of
02396                                    // the given argument are set that
02397                                    // make sense for a given space
02398                                    // dimension
02399   Assert ((refinement_case & RefinementPossibilities<dim>::isotropic_refinement) ==
02400           refinement_case,
02401           ExcInvalidRefinementCase (refinement_case));
02402 }
02403 
02404 
02405 
02406 template <int dim>
02407 inline
02408 RefinementCase<dim>::RefinementCase (const unsigned char refinement_case)
02409                 :
02410                 value (refinement_case)
02411 {
02412                                    // check that only those bits of
02413                                    // the given argument are set that
02414                                    // make sense for a given space
02415                                    // dimension
02416   Assert ((refinement_case & RefinementPossibilities<dim>::isotropic_refinement) ==
02417           refinement_case,
02418           ExcInvalidRefinementCase (refinement_case));
02419 }
02420 
02421 
02422 
02423 template <int dim>
02424 inline
02425 RefinementCase<dim>::operator unsigned char () const
02426 {
02427   return value;
02428 }
02429 
02430 
02431 
02432 template <int dim>
02433 inline
02434 RefinementCase<dim>
02435 RefinementCase<dim>::operator | (const RefinementCase<dim> &r) const
02436 {
02437   return RefinementCase<dim>(static_cast<unsigned char> (value | r.value));
02438 }
02439 
02440 
02441 
02442 template <int dim>
02443 inline
02444 RefinementCase<dim>
02445 RefinementCase<dim>::operator & (const RefinementCase<dim> &r) const
02446 {
02447   return RefinementCase<dim>(static_cast<unsigned char> (value & r.value));
02448 }
02449 
02450 
02451 
02452 template <int dim>
02453 inline
02454 RefinementCase<dim>
02455 RefinementCase<dim>::operator ~ () const
02456 {
02457   return RefinementCase<dim>(static_cast<unsigned char> (
02458                                (~value) & RefinementPossibilities<dim>::isotropic_refinement));
02459 }
02460 
02461 
02462 
02463 
02464 template <int dim>
02465 inline
02466 std::size_t
02467 RefinementCase<dim>::memory_consumption ()
02468 {
02469   return sizeof(RefinementCase<dim>);
02470 }
02471 
02472 
02473 
02474 template <int dim>
02475 template <class Archive>
02476 void RefinementCase<dim>::serialize (Archive &ar,
02477                                      const unsigned int)
02478 {
02479   // serialization can't deal with bitfields, so copy from/to a full sized
02480   // unsigned char
02481   unsigned char uchar_value = value;
02482   ar & uchar_value;
02483   value = uchar_value;
02484 }
02485 
02486 
02487 
02488 
02489 template <>
02490 inline
02491 Point<1>
02492 GeometryInfo<1>::unit_cell_vertex (const unsigned int vertex)
02493 {
02494   Assert (vertex < vertices_per_cell,
02495           ExcIndexRange (vertex, 0, vertices_per_cell));
02496 
02497   return Point<1>(static_cast<double>(vertex));
02498 }
02499 
02500 
02501 
02502 template <>
02503 inline
02504 Point<2>
02505 GeometryInfo<2>::unit_cell_vertex (const unsigned int vertex)
02506 {
02507   Assert (vertex < vertices_per_cell,
02508           ExcIndexRange (vertex, 0, vertices_per_cell));
02509 
02510   return Point<2>(vertex%2, vertex/2);
02511 }
02512 
02513 
02514 
02515 template <>
02516 inline
02517 Point<3>
02518 GeometryInfo<3>::unit_cell_vertex (const unsigned int vertex)
02519 {
02520   Assert (vertex < vertices_per_cell,
02521           ExcIndexRange (vertex, 0, vertices_per_cell));
02522 
02523   return Point<3>(vertex%2, vertex/2%2, vertex/4);
02524 }
02525 
02526 
02527 
02528 template <int dim>
02529 inline
02530 Point<dim>
02531 GeometryInfo<dim>::unit_cell_vertex (const unsigned int)
02532 {
02533   Assert(false, ExcNotImplemented());
02534 
02535   return Point<dim> ();
02536 }
02537 
02538 
02539 
02540 template <>
02541 inline
02542 unsigned int
02543 GeometryInfo<1>::child_cell_from_point (const Point<1> &p)
02544 {
02545   Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
02546 
02547   return (p[0] <= 0.5 ? 0 : 1);
02548 }
02549 
02550 
02551 
02552 template <>
02553 inline
02554 unsigned int
02555 GeometryInfo<2>::child_cell_from_point (const Point<2> &p)
02556 {
02557   Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
02558   Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
02559 
02560   return (p[0] <= 0.5 ?
02561           (p[1] <= 0.5 ? 0 : 2) :
02562           (p[1] <= 0.5 ? 1 : 3));
02563 }
02564 
02565 
02566 
02567 template <>
02568 inline
02569 unsigned int
02570 GeometryInfo<3>::child_cell_from_point (const Point<3> &p)
02571 {
02572   Assert ((p[0] >= 0) && (p[0] <= 1), ExcInvalidCoordinate(p[0]));
02573   Assert ((p[1] >= 0) && (p[1] <= 1), ExcInvalidCoordinate(p[1]));
02574   Assert ((p[2] >= 0) && (p[2] <= 1), ExcInvalidCoordinate(p[2]));
02575 
02576   return (p[0] <= 0.5 ?
02577           (p[1] <= 0.5 ?
02578            (p[2] <= 0.5 ? 0 : 4) :
02579            (p[2] <= 0.5 ? 2 : 6)) :
02580           (p[1] <= 0.5 ?
02581            (p[2] <= 0.5 ? 1 : 5) :
02582            (p[2] <= 0.5 ? 3 : 7)));
02583 }
02584 
02585 
02586 template <int dim>
02587 inline
02588 unsigned int
02589 GeometryInfo<dim>::child_cell_from_point (const Point<dim> &)
02590 {
02591   Assert(false, ExcNotImplemented());
02592 
02593   return 0;
02594 }
02595 
02596 
02597 
02598 template <>
02599 inline
02600 Point<1>
02601 GeometryInfo<1>::cell_to_child_coordinates (const Point<1>         &p,
02602                                             const unsigned int      child_index,
02603                                             const RefinementCase<1> refine_case)
02604 
02605 {
02606   Assert (child_index < 2,
02607           ExcIndexRange (child_index, 0, 2));
02608   Assert (refine_case==RefinementCase<1>::cut_x,
02609           ExcInternalError());
02610 
02611   return p*2.0-unit_cell_vertex(child_index);
02612 }
02613 
02614 
02615 
02616 template <>
02617 inline
02618 Point<2>
02619 GeometryInfo<2>::cell_to_child_coordinates (const Point<2>         &p,
02620                                             const unsigned int      child_index,
02621                                             const RefinementCase<2> refine_case)
02622 
02623 {
02624   Assert (child_index < GeometryInfo<2>::n_children(refine_case),
02625           ExcIndexRange (child_index, 0, GeometryInfo<2>::n_children(refine_case)));
02626 
02627   Point<2> point=p;
02628   switch (refine_case)
02629     {
02630       case RefinementCase<2>::cut_x:
02631             point[0]*=2.0;
02632             if (child_index==1)
02633               point[0]-=1.0;
02634             break;
02635       case RefinementCase<2>::cut_y:
02636             point[1]*=2.0;
02637             if (child_index==1)
02638               point[1]-=1.0;
02639             break;
02640       case RefinementCase<2>::cut_xy:
02641             point*=2.0;
02642             point-=unit_cell_vertex(child_index);
02643             break;
02644       default:
02645             Assert(false, ExcInternalError());
02646     }
02647 
02648   return point;
02649 }
02650 
02651 
02652 
02653 template <>
02654 inline
02655 Point<3>
02656 GeometryInfo<3>::cell_to_child_coordinates (const Point<3>         &p,
02657                                             const unsigned int      child_index,
02658                                             const RefinementCase<3> refine_case)
02659 
02660 {
02661   Assert (child_index < GeometryInfo<3>::n_children(refine_case),
02662           ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case)));
02663 
02664   Point<3> point=p;
02665                                    // there might be a cleverer way to do
02666                                    // this, but since this function is called
02667                                    // in very few places for initialization
02668                                    // purposes only, I don't care at the
02669                                    // moment
02670   switch (refine_case)
02671     {
02672       case RefinementCase<3>::cut_x:
02673             point[0]*=2.0;
02674             if (child_index==1)
02675               point[0]-=1.0;
02676             break;
02677       case RefinementCase<3>::cut_y:
02678             point[1]*=2.0;
02679             if (child_index==1)
02680               point[1]-=1.0;
02681             break;
02682       case RefinementCase<3>::cut_z:
02683             point[2]*=2.0;
02684             if (child_index==1)
02685               point[2]-=1.0;
02686             break;
02687       case RefinementCase<3>::cut_xy:
02688             point[0]*=2.0;
02689             point[1]*=2.0;
02690             if (child_index%2==1)
02691               point[0]-=1.0;
02692             if (child_index/2==1)
02693               point[1]-=1.0;
02694             break;
02695       case RefinementCase<3>::cut_xz:
02696                                              // careful, this is slightly
02697                                              // different from xy and yz due to
02698                                              // differnt internal numbering of
02699                                              // children!
02700             point[0]*=2.0;
02701             point[2]*=2.0;
02702             if (child_index/2==1)
02703               point[0]-=1.0;
02704             if (child_index%2==1)
02705               point[2]-=1.0;
02706             break;
02707       case RefinementCase<3>::cut_yz:
02708             point[1]*=2.0;
02709             point[2]*=2.0;
02710             if (child_index%2==1)
02711               point[1]-=1.0;
02712             if (child_index/2==1)
02713               point[2]-=1.0;
02714             break;
02715       case RefinementCase<3>::cut_xyz:
02716             point*=2.0;
02717             point-=unit_cell_vertex(child_index);
02718             break;
02719       default:
02720             Assert(false, ExcInternalError());
02721     }
02722 
02723   return point;
02724 }
02725 
02726 
02727 
02728 template <int dim>
02729 inline
02730 Point<dim>
02731 GeometryInfo<dim>::cell_to_child_coordinates (const Point<dim>         &/*p*/,
02732                                               const unsigned int        /*child_index*/,
02733                                               const RefinementCase<dim> /*refine_case*/)
02734 
02735 {
02736   AssertThrow (false, ExcNotImplemented());
02737   return Point<dim>();
02738 }
02739 
02740 
02741 
02742 template <>
02743 inline
02744 Point<1>
02745 GeometryInfo<1>::child_to_cell_coordinates (const Point<1>         &p,
02746                                             const unsigned int      child_index,
02747                                             const RefinementCase<1> refine_case)
02748 
02749 {
02750   Assert (child_index < 2,
02751           ExcIndexRange (child_index, 0, 2));
02752   Assert (refine_case==RefinementCase<1>::cut_x,
02753           ExcInternalError());
02754 
02755   return (p+unit_cell_vertex(child_index))*0.5;
02756 }
02757 
02758 
02759 
02760 template <>
02761 inline
02762 Point<3>
02763 GeometryInfo<3>::child_to_cell_coordinates (const Point<3>         &p,
02764                                             const unsigned int      child_index,
02765                                             const RefinementCase<3> refine_case)
02766 
02767 {
02768   Assert (child_index < GeometryInfo<3>::n_children(refine_case),
02769           ExcIndexRange (child_index, 0, GeometryInfo<3>::n_children(refine_case)));
02770 
02771   Point<3> point=p;
02772                                    // there might be a cleverer way to do
02773                                    // this, but since this function is called
02774                                    // in very few places for initialization
02775                                    // purposes only, I don't care at the
02776                                    // moment
02777   switch (refine_case)
02778     {
02779       case RefinementCase<3>::cut_x:
02780             if (child_index==1)
02781               point[0]+=1.0;
02782             point[0]*=0.5;
02783             break;
02784       case RefinementCase<3>::cut_y:
02785             if (child_index==1)
02786               point[1]+=1.0;
02787             point[1]*=0.5;
02788             break;
02789       case RefinementCase<3>::cut_z:
02790             if (child_index==1)
02791               point[2]+=1.0;
02792             point[2]*=0.5;
02793             break;
02794       case RefinementCase<3>::cut_xy:
02795             if (child_index%2==1)
02796               point[0]+=1.0;
02797             if (child_index/2==1)
02798               point[1]+=1.0;
02799             point[0]*=0.5;
02800             point[1]*=0.5;
02801             break;
02802       case RefinementCase<3>::cut_xz:
02803                                              // careful, this is slightly
02804                                              // different from xy and yz due to
02805                                              // differnt internal numbering of
02806                                              // children!
02807             if (child_index/2==1)
02808               point[0]+=1.0;
02809             if (child_index%2==1)
02810               point[2]+=1.0;
02811             point[0]*=0.5;
02812             point[2]*=0.5;
02813             break;
02814       case RefinementCase<3>::cut_yz:
02815             if (child_index%2==1)
02816               point[1]+=1.0;
02817             if (child_index/2==1)
02818               point[2]+=1.0;
02819             point[1]*=0.5;
02820             point[2]*=0.5;
02821             break;
02822       case RefinementCase<3>::cut_xyz:
02823             point+=unit_cell_vertex(child_index);
02824             point*=0.5;
02825             break;
02826       default:
02827             Assert(false, ExcInternalError());
02828     }
02829 
02830   return point;
02831 }
02832 
02833 
02834 
02835 template <>
02836 inline
02837 Point<2>
02838 GeometryInfo<2>::child_to_cell_coordinates (const Point<2>         &p,
02839                                             const unsigned int      child_index,
02840                                             const RefinementCase<2> refine_case)
02841 {
02842   Assert (child_index < GeometryInfo<2>::n_children(refine_case),
02843           ExcIndexRange (child_index, 0, GeometryInfo<2>::n_children(refine_case)));
02844 
02845   Point<2> point=p;
02846   switch (refine_case)
02847     {
02848       case RefinementCase<2>::cut_x:
02849             if (child_index==1)
02850               point[0]+=1.0;
02851             point[0]*=0.5;
02852             break;
02853       case RefinementCase<2>::cut_y:
02854             if (child_index==1)
02855               point[1]+=1.0;
02856             point[1]*=0.5;
02857             break;
02858       case RefinementCase<2>::cut_xy:
02859             point+=unit_cell_vertex(child_index);
02860             point*=0.5;
02861             break;
02862       default:
02863             Assert(false, ExcInternalError());
02864     }
02865 
02866   return point;
02867 }
02868 
02869 
02870 
02871 template <int dim>
02872 inline
02873 Point<dim>
02874 GeometryInfo<dim>::child_to_cell_coordinates (const Point<dim>         &/*p*/,
02875                                               const unsigned int        /*child_index*/,
02876                                               const RefinementCase<dim> /*refine_case*/)
02877 {
02878   AssertThrow (false, ExcNotImplemented());
02879   return Point<dim>();
02880 }
02881 
02882 
02883 
02884 template <>
02885 inline
02886 bool
02887 GeometryInfo<1>::is_inside_unit_cell (const Point<1> &p)
02888 {
02889   return (p[0] >= 0.) && (p[0] <= 1.);
02890 }
02891 
02892 
02893 
02894 template <>
02895 inline
02896 bool
02897 GeometryInfo<2>::is_inside_unit_cell (const Point<2> &p)
02898 {
02899   return (p[0] >= 0.) && (p[0] <= 1.) &&
02900          (p[1] >= 0.) && (p[1] <= 1.);
02901 }
02902 
02903 
02904 
02905 template <>
02906 inline
02907 bool
02908 GeometryInfo<3>::is_inside_unit_cell (const Point<3> &p)
02909 {
02910   return (p[0] >= 0.) && (p[0] <= 1.) &&
02911          (p[1] >= 0.) && (p[1] <= 1.) &&
02912          (p[2] >= 0.) && (p[2] <= 1.);
02913 }
02914 
02915 template <>
02916 inline
02917 bool
02918 GeometryInfo<1>::is_inside_unit_cell (const Point<1> &p,
02919                                       const double eps)
02920 {
02921   return (p[0] >= -eps) && (p[0] <= 1.+eps);
02922 }
02923 
02924 
02925 
02926 template <>
02927 inline
02928 bool
02929 GeometryInfo<2>::is_inside_unit_cell (const Point<2> &p,
02930                                       const double eps)
02931 {
02932   const double l = -eps, u = 1+eps;
02933   return (p[0] >= l) && (p[0] <= u) &&
02934          (p[1] >= l) && (p[1] <= u);
02935 }
02936 
02937 
02938 
02939 template <>
02940 inline
02941 bool
02942 GeometryInfo<3>::is_inside_unit_cell (const Point<3> &p,
02943                                       const double eps)
02944 {
02945   const double l = -eps, u = 1.0+eps;
02946   return (p[0] >= l) && (p[0] <= u) &&
02947          (p[1] >= l) && (p[1] <= u) &&
02948          (p[2] >= l) && (p[2] <= u);
02949 }
02950 
02951 
02952 #endif // DOXYGEN
02953 DEAL_II_NAMESPACE_CLOSE
02954 
02955 #endif
02956 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Thu May 17 2012 20:04:34 by doxygen 1.7.3