00001
00002
00003
00004
00005
00006
00007
00008
00009
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 }
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
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
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 }
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
02396
02397
02398
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
02413
02414
02415
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
02480
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
02666
02667
02668
02669
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
02697
02698
02699
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> &,
02732 const unsigned int ,
02733 const RefinementCase<dim> )
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
02773
02774
02775
02776
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
02804
02805
02806
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> &,
02875 const unsigned int ,
02876 const RefinementCase<dim> )
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